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We present the first generalization of Navier-Stokes theory to relativity that satisfies all of the following 
properties: (a) the system coupled to Einstein’s equations is causal and strongly hyperbolic; (b) equilibrium 
states are stable; (c) all leading dissipative contributions are present, i.e., shear viscosity, bulk viscosity, and 
thermal conductivity; (d) nonzero baryon number is included; (e) entropy production is non-negative in the 
regime of validity of the theory; (f) all of the above hold in the nonlinear regime without any simplifying 
symmetry assumptions. These properties are accomplished using a generalization of Eckart’s theory 
containing only the hydrodynamic variables, so that no new extended degrees of freedom are needed as in 
Miiller-Israel-Stewart theories. Property (b), in particular, follows from a more general result that we also 
establish, namely, sufficient conditions that when added to stability in the fluid’s rest frame imply stability 
in any reference frame obtained via a Lorentz transformation All of our results are mathematically 
rigorously established. The framework presented here provides the starting point for systematic 
investigations of general-relativistic viscous phenomena in neutron star mergers. 
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I. INTRODUCTION 


Relativistic fluid dynamics has been successfully used as 
an effective description of long wavelength, long time 
phenomena in a multitude of different physical systems, 
ranging from cosmology [1] to astrophysics [2] and also 
high-energy nuclear physics [3]. In the latter, relativistic 
viscous fluid dynamics has played an essential role in the 
description of the dynamical evolution of the quark-gluon 
plasma formed in ultrarelativistic heavy-ion collisions [4] 
and also in the quantitative extraction of its transport 
properties (see, for instance, Ref. [5]). More recently, with 
the observation of binary neutron star mergers [6-8], the 
modeling of the different dynamical stages experienced by 
the hot and dense matter formed in these collisions requires 
extending of our current understanding of viscous fluids 
toward the strong gravity regime where general relativistic 
effects are important (see, e.g., Refs. [9-14]. 
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The ubiquitousness of fluid dynamics stems from the 
existence of general conservation laws (such as energy, 
momentum, and baryon number) and their consequences to 
systems where there is a large separation of scales, such that 
the macroscopic behavior of conserved quantities can be 
understood without precise knowledge of all the details that 
govern the system’s underlying microscopic properties 
[15]. Ideal fluid dynamics is the extreme situation where 
dissipative effects are neglected and the theory’s basic 
properties in this limit are reasonably well understood, both 
in a fixed background as well as when coupling to 
Einstein’s equations is taken into account [2,16,17]. We 
remark that because all sources of dissipation relevant for 
our discussion stem from bulk viscosity, shear viscosity, 
and heat conduction, and following standard practice in the 
field [3], we will use the terms viscous fluid and dissipative 
fluid interchangeably. In particular, other sources of dis- 
sipation, such as anomalous dissipation [18,19], will not be 
discussed. 

When dissipative effects are taken into account, the 
behavior of fluids is far less understood (unless stated 
otherwise, fluids, hydrodynamics, etc., henceforth mean 
relativistic fluids, relativistic hydrodynamics, etc.), despite 
the importance of viscous dissipation in cutting-edge scien- 
tific experiments such as in studies of the quark-gluon plasma 
or their expected relevance for neutron star mergers, as 
mentioned above. Historically, a stumbling block has been 
the difficulty of modeling dissipative phenomena while 
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preserving causality. Causality is a central postulate in 
special and general relativity, stating that the speed which 
information can propagate in any system cannot be larger 
than the speed of light [20]. This implies that a solution to the 
equations of motion at a given spacetime point x is 
completely determined by the spacetime region that is in 
the past of and causally connected to x [20-22]. Of course, 
this property must hold in relativity regardless of whether 
dissipation is present or not [20]. While causality is typically 
not an issue for most matter models under reasonable 
assumptions [21], including the case of ideal fluids [23], 
ensuring causality of fluid theories in the presence of 
dissipation turned out to be a major challenge [2]. 

The challenges one encounters when modeling fluids 
with dissipation, however, are not restricted to enforcing 
causality. Another hallmark property of dissipative fluids is 
stability. By this we mean that perturbations of a system 
that is in thermodynamic equilibrium should decay in time. 
This expresses the basic intuition that if dissipation is 
present, the system will dissipate energy and, consequently, 
small deviations from equilibrium will be damped, leading 
the dynamics to return to equilibrium within some char- 
acteristic timescale. Naturally, in order to implement this 
idea in a given formalism one needs to specify what is 
meant by equilibrium and perturbations. We will consider 
homogeneous (nonrotating) equilibrium states and our 
perturbations will be plane-wave solutions to the equations 
of motion linearized about such homogeneous states. 
Although this is not the most general definition of equi- 
librium [24], it captures the most basic intuition about how 
deviations from equilibrium should behave in a dissipative 
theory and, consequently, in practice this has been the 
definition most often used in the literature [25,26]. Like 
causality, stability is a property that is difficult to incor- 
porate in theories of relativistic fluids with dissipation. 

Aside from causality and stability, a third fundamental 
property required for a theory of relativistic viscous fluids 
is that the equations of motion be locally well posed. This 
means that given initial conditions, there must exist one and 
only one solution to the equations of motion taking the 
prescribed initial conditions [27] and defined for some time 
[28]. Physically, this means that the system has a well- 
defined evolution determined by the initial conditions. Like 
causality, local well posedness is a property required of any 
field theory [2,20—22], but we emphasize it here since, also 
like causality, this is a property that is difficult to achieve in 
theories of fluids with dissipation. 

Needless to say, there is little use for a theory of fluids 
that is causal, stable, and locally well posed if it is not able 
to make connections with real physical phenomena. Thus, a 
theory of relativistic viscous fluids must in addition be 
suitable for empirical studies. This means, at the least, that 
the theory must agree with well-established physical facts, 
but also that one needs to be able to extract quantitative 
predictions from such a theory. 
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The interplay between theory and experiment is, of 
course, at the heart of physics. In the context of relativistic 
fluid dynamics, such interplay has been heavily guided by 
complex numerical simulations [3]. Moreover, it is clear 
that simulations will continue to be at the center of 
developments in the field, particularly when it comes to 
the investigations of viscous effects in neutron star mergers. 
In this regard, while there is no one-size-fits-all approach 
for implementing numerical simulations of general relativ- 
istic systems [2,29], in the numerical general relativity 
community one concept that has been very important for 
the construction of numerical algorithms is that of strong 
hyperbolicity [30]. This means that the principal part of the 
equations of motion can be diagonalized; see Sec. V for 
details. Although a discussion of the role of strong hyper- 
bolicity in general relativistic numerical simulations is 
beyond the scope of this work (the reader can consult 
the above references for details), we stress that strong 
hyperbolicity is a highly desirable feature for numerical 
studies of general relativistic systems (see also Ref. [31] for 
more discussion on potential caveats of numerical 
simulations). 

In summary, a physically meaningful theory of relativ- 
istic viscous fluids must be (I) causal, (II) stable, and (II) 
locally well posed. In addition, it is highly desirable to have 
a theory that is ([V) strong hyperbolic. 

While property II is, by definition, concerned with the 
equations linearized about equilibrium in Minkowski back- 
ground, we emphasize that whenever referring to causality, 
local well posedness, and strong hyperbolicity, i.e., proper- 
ties I, II, and IV, we are always talking about the equations 
of motion in the full nonlinear regime. It is important to 
stress this point because a substantial body of theoretical 
work in relativistic viscous fluids is restricted to analyzing 
the equations linearized about equilibrium and, thus, the 
corresponding claims about causality, local well posedness, 
etc., are restricted to this particular, linearized-about-equi- 
librium case (see Sec. II B). Furthermore, for applications 
in general relativity (in particular the study of viscous 
effects in neutron star mergers), one is interested in the case 
where properties I-IV hold with dynamical coupling to 
Eintein’s equations (again, with exception of property ID). 

At this point, we should stress that when we say that a 
theory is causal, stable, etc., we do not mean it uncondi- 
tionally, but rather under a specific set of assumptions. 
Obviously, one is interested in cases where the assumptions 
are physically reasonable, even if they do not cover all 
cases of physical interest. For simplicity, however, in the 
remaining of this Introduction and in Sec. II, we avoid 
discussion of specific hypotheses. Thus, when we say that a 
certain theory is causal, stable, etc., we mean “causal under 
a specific set of assumptions,” and unless stated otherwise, 
it will be implicitly understood that the assumptions in 
question are of physical interest. An exception to this will 
be made only later in Sec. IIB, when we summarize the 
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extent to which different theories of viscous fluids satisfy 
one or more of the properties I-IV, since in this case 
mentioning the assumptions under which such theories 
fulfill some of these requirements will be important for 
comparison among them and also with our results. Even in 
this case, however, we will refer to those assumptions only 
at a high level (e.g., we will say that a certain property holds 
for nonzero shear viscosity but without specifying the 
precise range of nonzero values that is in fact required for 
the result to hold). We believe that this will suffice to give 
the reader a panoramic view of the state of affairs in the 
field. All the precise assumptions for the results that will be 
discussed can be found in the references we provide or, in 
the case of the results of this paper, in the remainder of 
the text. 

The goal of this work is to provide the first example of a 
theory of relativistic viscous fluids that simultaneously 
satisfies all the properties I-IV. All our results are math- 
ematically rigorous, hold with or without dynamical 
coupling to Einstein’s equations, are valid in the full 
nonlinear regime, and do not make any symmetry or 
simplifying assumptions. We establish these results without 
the need for additional (extended) variables (see Sec. IIB 
for details). 

Section II provides a more or less self-contained expo- 
sition of our results and how they fit within studies of 
relativistic fluids with viscosity. We hope that such an 
exposition will be helpful to readers interested in the 
subject here investigated but who are not necessarily 
specialists in all the topics covered by our methods. In 
order to keep our account as simple as possible, we carry 
out the discussion in Sec. II at a high level, writing few 
formulas and omitting several details, but we provide full 
references for interested readers. More precisely, in Sec. II 
A, we discuss some important concepts underlying the 
investigation of relativistic viscous fluids. None of the ideas 
discussed in Sec. II A are new, but they play a key role in 
our constructions. Therefore, it is convenient to revisit such 
ideas here. In Sec. IIB, we review the state of affairs in the 
field regarding properties I-IV. This review is not intended 
to be exhaustive; rather, our goal is to provide enough 
context for our results. Finally, in Sec. IIC, we provide a 
summary and discussion of our results. Specialists might 
skip Sec. II without compromising understanding (although 
some specialists might still be interested in some aspects of 
the discussion in Sec. II C). 

Definitions.—The spacetime metric g,,, has a mostly plus 
signature (— + ++). Greek indices run from 0 to 3, latin 
indices from | to 3. The spacetime covariant derivative is 
denoted as V,,. We use natural units, c = h = kg = 1. 


A. Organization of the paper 


This paper is organized as follows. In Sec. II we provide 
an overview of our results and the context surrounding 
them. In Sec. II, we formulate a generalization of Navier- 
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Stokes (NS) theory using the Bemfica-Disconzi-Noronha- 
Kovtun (BDNK) formalism [32—35]. In Sec. IV, we provide 
necessary and sufficient conditions that must be fulfilled by 
the parameters of the theory for causality to hold. In Sec. V, 
we prove that the full nonlinear system of equations in 
general relativity is strongly hyperbolic, the solutions are 
unique, and the initial-value problem is well posed in 
general relativity. A new theorem concerning the linear 
stability properties of relativistic fluids in flat spacetime is 
given in Sec. VI. We employ this theorem in Sec. VII to 
obtain conditions that ensure that the new theory presented 
here is stable. The rigorous mathematical proofs of 
Theorem I, Proposition I, Theorem II, and Theorem II 
are found in the Appendixes A, B, C, and D, respectively. 
Our conclusions and outlook can be found in Sec. VII. 


Il. BACKGROUND AND DISCUSSION 


A. Definition of out-of-equilibrium variables: 
Hydrodynamic frames 


In the modern perspective, relativistic fluid dynamics is 
understood as an effective theory for the evolution of 
conserved densities, such as the energy-momentum tensor 
T’”. (We could include, in this introductory part, other 
conserved quantities such as the baryon current J“ and 
those associated with higher moments. In fact, conservation 
of J” will be implicitly understood later in the discussion of 
Secs. II A and ITB and thereafter since we will often refer to 
the presence of a chemical potential. For simplicity, 
however, we will often refer only to 7” in this part, since 
this will suffice for the aspects we want to highlight.) To say 
that 7“” is conserved means that 


V,T" = 0, 


which provides equations of motion governing the dynam- 
ics of the fluid. 

The energy-momentum tensor 7’” is understood as the 
expectation value of the microscopic quantum operator 7”, 
which is an observable that can be defined for any non- 
equilibrium state. In equilibrium, the state of the system can 
be parametrized by the temperature T,,, the flow velocity 
ue (observe that this is the four-velocity of the fluid, 
although we will often refer to it simply as the velocity; the 
fluid velocity is always assumed to be normalized; see 
Sec. IIB), and the chemical potential yg. One of the 
assumptions that forms the basis of a fluid dynamics 
description is that for states not very far from equilibrium, 
the physical observable T#” = (7*”) can still be para- 
metrized in terms of a “temperature” 7, a “flow velocity” 
u“, and a “chemical potential” yw that reduce to i tea and 
Meg in equilibrium. We write quotation marks to emphasize 
the fact that the quantities T, u“, and yw have no first- 
principles microscopic definitions. Therefore, while it is 
useful to interpret T, uw“, and mw as out-of-equilibrium 
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macroscopic temperature, velocity, and chemical potential, 
since they are close to Tg, Ueg, and Meg and reduce to the 
latter in equilibrium, we should ultimately understand T, 
u“, and yas auxiliary variables that are used to parametrize 
the physical observable T“”—the latter enjoying a first- 
principles, microscopic definition even when the system is 
out of equilibrium. 

It follows that there exists an ambiguity in the definition 
of the out-of-equilibrium quantities 7, uv“, and yu, since there 
are different ways of parametrizing T’” subject to the 
constraint that one recovers the unambiguous parametriza- 
tion in terms of Tg, Ueg, and Meg in equilibrium. In other 
words, different out-of-equilibrium choices of T, u”, and hu 
to parametrize T’” are allowed as long as they agree in 
equilibrium. This is sometimes expressed by saying that T, 
u“, and yw correspond to a “fictitious” temperature, flow 
velocity, and chemical potential [2,36]. 

A particular choice of parametrization of T’” in terms of 
T, u“, and yw has been historically called a choice of a 
hydrodynamic frame, or simply frame. (It is a bit unfortu- 
nate the word “frame” has also other meanings in relativity 
theory, e.g., reference frames related by a Lorentz trans- 
formation, frames in a tetrad formalism, null frames, or a 
local rest frame (LRF), etc. However, all these different 
meanings can be distinguished from the context.) A choice 
of frame is, therefore, a definition of what one means by 
temperature, velocity, and chemical potential out of equi- 
librium. Consequently, a choice of frame is always involved 
whenever we describe a fluid out of equilibrium in terms of 
temperature, velocity, and chemical potential. This is still 
the case even if further, extended variables are introduced 
(see Sec. IIB for the notion of extended variables). The 
notion of hydrodynamic frame and how it represents a 
choice of out-of-equilibrium variables is discussed exten- 
sively in the literature. An incomplete list is given by 
Refs. [34,36—-51]. References [34,48], in particular, contain 
a detailed discussion of the topic. 

Observe that once T”” is cast in terms of 7, uv”, and p, the 
energy-momentum conservation equations V,,7"” = 0 can 
be equivalently written as evolution equations for those 
quantities. We also remark that one can choose other 
thermodynamic quantities, e.g., the energy density or the 
pressure, to parametrize T””, and we will in fact do so later on 
in the paper. Of course, not all thermodynamic scalars are 
independent; they are connected by the first law of thermo- 
dynamics and a prescription of an equation of state [2]. 
Obviously, the nonuniqueness in the definition of the 
variables used to parametrize T’” out of equilibrium remains 
if we choose a parametrization in terms of other thermody- 
namic variables such as the energy density, pressure, etc. 

In order to pass from this qualitative argument about the 
ambiguity of T, uv“, and » away from equilibrium to a more 
precise assessment of such ambiguity, one needs to be more 
specific about how one formalizes the idea that fluid 
dynamics arises as a long time, long wavelength limit of 
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an underlying microscopic theory, i.e., as a description of 
the macroscopic dynamics of the system for small devia- 
tions from equilibrium. Such a formalization can be 
accomplished in the framework of the so-called gradient 
expansion, which was used a century ago by Chapman and 
Enskog in the derivation of fluid dynamics from the 
(nonrelativistic) Boltzmann equation and that has since 
then been adapted to the relativistic setting [38]. We remark 
that the gradient expansion is not the only way to formalize 
the idea that fluid dynamics is an effective description that 
emerges from a more fundamental microscopic behavior; 
see Sec. IIB for a discussion of ideas involving the so- 
called moment expansion and holographic techniques. 
Nevertheless, the gradient expansion, while not fundamen- 
tal, is a very convenient and powerful formalism based on 
effective field theory ideas that allows one to track how 
different parametrizations of T“” lead to different fluid 
descriptions. 

The gradient expansion is based on the idea that one can 
write 


T” = O(1)+ O(0) + O(07) + ---, 





where O(0") denotes terms with n derivatives of T, wu“, and 
y [so, e.g., O(0”) involves both terms of the form 077 and 
OT Op, etc.] and O(1) corresponds to the terms that reduce 
to Téq, the energy-momentum tensor parametrized in terms 
Of Tas Ueg, and Meg: Schematically, this is an expansion in 
powers of the Knudsen number Kn ~ @pic,90, i-e., the ratio 
between the relevant microscopic scale (pico and the 
inverse macroscopic scale L, associated with the derivative 
of the hydrodynamic fields. In this sense, the gradient 
expansion corresponds to the well-known Knudsen number 
expansion used in the description of kinetic systems 
[38,39]. In particular, since the expansion truncated at 
O(1) corresponds to ideal hydrodynamics, viscous con- 
tributions require considering at least O(O) terms, which is 
consistent with the basic intuition that dissipation is a 
phenomenon associated with deviations from equilibrium. 

In order to construct a fluid theory out of the gradient 
expansion, one truncates it at a certain order. This trunca- 
tion necessarily defines a scale at which the effective 
description is supposed to be valid, with higher-order 
effects encoded by the terms neglected in the expansion 
which are considered outside the limit of validity of the 
truncated theory. Aside from the truncation order, one also 
needs to specify the constitutive relations, i.e., the specific 
form of each term O(0") in terms of T,u“,, up to the 
truncation order (see Secs. IIB and III for examples). By 
specifying the truncation order and the constitutive rela- 
tions, one is in fact defining what is meant by T, vu", and yu 
out of equilibrium; i.e., one is making a choice of hydro- 
dynamic frame. 

Different frame choices, therefore, correspond to differ- 
ent effective descriptions of the same truncated theory. At 
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this point, it seems almost unnecessary to talk about frames, 
and one might be tempted to simply say that one has 
distinct theories of fluids. The key word here, however, is 
effective. Indeed, when we consider two distinct constitu- 
tive relations truncated at a given order, 

T! =T!(T,u*,n) and TY = T!"(P, a, 7), 
one obviously has different fluid theories: the equations of 
motion V,7’” =O and Va" = (0 are not the same. 
Consequently (upon writing these conservation laws in 
terms of T, u%, ~ and T, au, ji, respectively), the quantities 
T,u*,u and T, it”, ji satisfy different evolution equations 
and, thus, cannot represent the same definition of temper- 
ature, fluid velocity, and chemical potential. However, 
one needs to keep in mind that the temperature, flow 
velocity, and chemical potential are not fundamental 
quantities, whereas the energy-momentum tensor is (it 
does have a first-principles definition). Thus, T’”(T, u%, 1) 
and T"’(T, a", ji) differ because they represent distinct 
coarse-grained or low-energy limits of the actual, micro- 
scopically uniquely defined, energy-momentum tensor. 
Therefore, the language of frames signals the key fact that 
one is always considering one possible effective description 
among many. 

Summarizing, there exists an intrinsic ambiguity in how 
one parametrizes T“” in terms of out-of-equilibrium tem- 
perature 7, velocity u, and chemical potential w. Such 
ambiguity simply expresses the fact that these quantities do 
not have first-principles microscopic definitions away from 
equilibrium. What is not ambiguous away from equilibrium 
is the definition of T“”. One resolves this ambiguity by 
choosing a definition of T, wu", w. Such a choice is known as 
a frame choice. Different parametrizations of T””, therefore, 
correspond to different frame choices. Not all frame 
choices, however, are equally useful. In our work, we 
explore suitable definitions of temperature, flow velocity, 
and chemical potential to construct effective theories 
describing fluids that lead to sensible theories in terms 
of satisfying properties I-IV. 

At this point, the attentive reader will probably have 
noticed that much of the above discussion does not depend 
on relativistic principles. In other words, the fact that there 
is no first-principles definition of out-of-equilibrium quan- 
tities such as temperature, flow velocity, and chemical 
potential applies to nonrelativistic theories as well. In the 
nonrelativistic setting, however, there exists a highly 
successful theory of dissipative (Newtonian) fluids, 
namely, the Navier-Stokes-Fourier theory. In light of its 
success, it is fair to say that for all practical purposes, one 
can take the definitions of out-of-equilibrium quantities in 
the Navier-Stokes-Fourier theory as the correct ones in a 
nonrelativistic context. Had an equivalently successful 
theory of relativistic viscous fluids been available (where 
success would in particular incorporate properties I-IV), 
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we could similarly take the definitions of out-of-equilib- 
rium quantities in such a theory as the correct ones for all 
practical purposes. Nevertheless, as we explain in the next 
section, there is not, at the moment, a theory of relativistic 
viscous fluids that can claim such a level of success. Hence, 
exploring how different frame choices can lead to different 
fluid descriptions becomes a topic of uttermost interest (see 
Sec. ITC). 


B. Brief overview of viscous theories 


The first proposal for a relativistic viscous fluid theory 
was done by Eckart [52] in 1940, with a closely related 
formulation by Landau and Lifshitz [15] in the 1950s. In 
these works, the authors postulated a form for the energy- 
momentum tensor (and also of the baryon current J“, but, 
as in the previous section, here we simplify the discussion 
by focusing on T“” only) based on ideas from thermody- 
namics and following a covariant generalization of the 
nonrelativistic Navier-Stokes-Fourier theory. For example, 
in Eckart’s theory, one has 


Tekan = €UMUY + (P —CV u*) AM” + gh’ + q’u" —2no™, 


where e, T, and uw“ are the (out-of-equilibrium) energy 
density [53], temperature, and velocity of the fluid, with the 
latter normalized [54] by uu, = —1, A” = g” + u*u’ is 
the projection onto the space orthogonal to wu“, P is the 
equilibrium pressure (see below) given by an equation of 
state (the choice of which depends on the nature of the 
fluid; for example, for a conformal fluid one has P = 5), € 
is the coefficient of bulk viscosity, 7 is the coefficient of 
shear viscosity, g, = —KT(AjV,InT + u’V,u,) repre- 
sents energy diffusion, with x being the coefficient of heat 
conduction, and of” = A“ Vall is the shear tensor, with 
ANp = 5 (AGA, + AAG —FAMA,y) (so Ate projects a 
two-tensor on the space of two-tensors traceless and 
orthogonal to wv“). In the absence of viscous effects, when 
€ =n =x = 0, one recovers the energy-momentum tensor 
of an ideal fluid. 

According to the standard physical interpretation of the 
energy-momentum tensor of a fluid, the fluid’s total 
pressure is given by 5A,,T™. It is convenient to write 
the total pressure as a sum of an “equilibrium” part, which 
is assumed to be given by an equation of state whose 
functional form follows that assigned to the fluid in the 
limit when viscous effects are absent, and a “nonequili- 
brium” part that contains explicitly the viscous contribu- 
tions. In the case of Ty... the latter is given by —CV_,u". 
This term clearly illustrates the fact that only terms of first 
order in Knudsen number were kept in this case because 
¢/P gives the relevant microscopic length scale associated 
with particle-number changing processes, while V,,u/ 
accounts for the inverse length scale associated with the 
gradient of the hydrodynamic fields. 
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As said, Eckart and Landau and Lifshitz were seeking a 
covariant version of the nonrelativistic Navier-Stokes 
equation compatible with thermodynamic principles, most 
notably, the second law of thermodynamics; i.e., their 
choice of 7“” ensured that entropy production (for a 
suitable definition of out-of-equilibrium entropy) is non- 
negative. From a modern perspective, however, these 
theories are better understood as effective theories that 
arise from a gradient expansion truncated at first order and 
with a specific choice of hydrodynamic frame, i.e., a 
specific choice of constitutive relation that parametrizes 
the energy-momentum tensor in terms of out-of-equilib- 
rium variables. In fact, it is possible to show that the Eckart 
and Landau and Lifshitz theories can be obtained from 
kinetic theory as an expansion in gradients truncated at first 
order [38]. Constraints on the coefficients that appear in 
such truncated series are found by imposing the second law 
of thermodynamics. In accordance with the notion of 
hydrodynamic frames, the specific choices that lead to 
the theories of Eckart and of Landau and Lifshitz are known 
in the literature as the Eckart and Landau and Lifshitz 
frames [2]. One can immediately see that other frame 
choices are possible for an energy-momentum tensor 
truncated at first order upon noticing that Tojo, does 
not contain all possible terms that are linear in derivatives of 
T, u“, and ~—terms that are allowed in a truncation at first 
order. Theories arising from a gradient expansion truncated 
at first order are known as first-order theories. The Eckart 
and Landau theories are, thus, examples of first-order 
theories. 

The Eckart and Landau and Lifshitz theories are very 
intuitive and natural at first sight. They correspond to 
immediate covariant generalizations of the nonrelativistic 
Navier-Stokes-Fourier theory (in fact, they recover it in the 
nonrelativistic limit), satisfy the second law of thermody- 
namics, preserve many features present in the ideal case 
(e.g., the energy density is recovered from the energy- 
momentum tensor by double contraction with the velocity), 
are relatively simple, and, as already said, can be derived 
from kinetic theory. Yet, they are remarkably at odds with 
fundamental physical principles in that they are known to 
violate causality and are unstable [25,55]. Consequently, 
the Eckart and Landau and Lifshitz theories cannot be taken 
as viable theories of relativistic viscous fluids. In fact, a 
large class of first-order theories, of which the theories of 
Eckart and of Landau and Lifshitz are particular cases, are 
known to be acausal and unstable [25]. One naturally 
wonders what are the root causes of the failures of these 
theories, especially when at first sight they look very 
intuitive. We return to this point in Sec. ITC. 

A different approach for the construction of relativistic 
viscous fluid theories was taken by Israel and Stewart in a 
series of works [36,37,56—58], adapting ideas developed by 
Miiller in the nonrelativistic setting [59]. The resulting 
theory is referred to Israel-Stewart or Miiller-Israel-Stewart 
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(MIS) theory, or sometimes simply Israel-Stewart theory. 
In the MIS theory, the energy momentum takes the form 


Tris = eu“u’ + (P+ TD) A” + OfuY + O’uk + ah”. 








The quantities I], z“”, and QO” represent the bulk viscosity, 
shear viscosity, and energy diffusion of the fluid, and are 
referred to as viscous fluxes. We see that T§.,,, corre- 
sponds to the choices where the bulk scalar H = —CV,,u", 
the shear-stress tensor is given by a¥” = —2no"”, and the 
energy diffusion reads Q! = q! = —xT(A;V,InT + 
u’V,u,)- In the MIS theory, however, the viscous fluxes 
are taken to be new variables on the same par as the 
“ordinary” variables T, u“, etc., (see below). Because 
II, 2”, QO” add to the number of variables, hence extending 
the state space, they are known as extended (thermody- 
namic) variables and theories that investigate extended 
variables are referred to as extended (thermodynamic) 
theories [60,61]. An important point to make (already 
alluded to earlier) is that one cannot dispense with a choice 
of hydrodynamic frame even in extended theories, since 
one still needs to make a definition of out-of-equilibrium 
temperature, flow velocity, and chemical potential. 

At this point, it is convenient to make the following 
definition. The variables 7, u“, » and those derived from 
them via the first law of thermodynamics and a choice of 
equation of state are known as hydrodynamic variables or 
fields. In other words, the hydrodynamic variables are the 
“ordinary” fields already present in the case of an ideal 
fluid (although the physical interpretation of these variables 
is not precisely the same as in the ideal fluid case; as 
discussed, the meaning of, e.g., temperature is different in 
or out of equilibrium). In this language, we can say that the 
Eckart and Landau and Lifshitz theories involve only the 
hydrodynamic variables, whereas the MIS theory involves 
both hydrodynamic and extended fields. In addition, the 
gradient expansion is always an expansion in the hydro- 
dynamic variables [62]. 

Because the MIS formalism introduces new variables in 
addition to the hydrodynamic fields, it also requires new 
equations of motion besides the standard conservation laws 
such as V,,Ty;3 = 0. The desired equations are postulated 
to be relaxation-type equations whose precise form is 
chosen so that entropy production is non-negative—where 
the entropy current is also extended from its usual form 
used in ideal fluids to include the extended variables 
II, z“”, O". For example, IT satisfies 


1 
rut V,T +1 = ~¢V,u" — eT, (= u) 


where ty is a relaxation time. See, e.g., Ref. [2] for the full 
set of equations satisfied by I], 2”, QO”, the form of the 
entropy current including these fields, and the derivation of 
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the equations of motion from the second law of thermo- 
dynamics [63]. 

The MIS theory enjoys the following good properties: 
the equations of motion are stable, thus satisfying property 
II, and their linearization about equilibrium states is causal, 
thus satisfying property I [24,64]. Also, it can, in certain 
limits, be derived from kinetic theory [2,57,58]. 

We next discuss three other theories of great interest that 
employ extended variables: Denicol-Niemi-Molnar-Rischke 
(DNMR), resumed  Baier-Romatschke-Son-Starinets- 
Stephanov (rBRSSS), and anisotropic hydrodynamics 
(AHYDRO) theories. The DNMR theory is an effective 
theory derived from kinetic theory via an expansion in 
moments [65]. The moment expansion goes back to Grad in 
his work on nonrelativistic fluids [66,67]. Applying this 
formalism to the relativistic Boltzmann equations, together 
with a new power-counting scheme involving Knudsen and 
inverse Reynolds number expansions, DNMR arrived at a 
set of equations for the hydrodynamic fields and a set of 
extended variables I, z“”, and Q* that represent the bulk 
viscosity, shear viscosity, and energy diffusion, similarly to 
the MIS equations. Also similar to the MIS equation is the 
fact that the equations satisfied by the viscous fluxes in the 
DNMR theory are relaxation-type equations. Despite their 
similarities, it is important to stress that the MIS and 
DNMR equations are not the same. 

The DNMR theory enjoys many good properties. It is 
stable and its linearization about equilibrium states [68] is 
causal [26]. When only bulk viscosity is present, the 
DNMR theory is causal, locally well posed, and strongly 
hyperbolic; these properties hold with and without dynami- 
cal coupling to Einstein’s equations [71]. When all viscous 
fluxes are present, but chemical potential is absent, the 
DNMR equations have recently been shown to be causal 
(again, with or without coupling to Einstein’s equations) 
[72] (see Refs. [26,73,74] for related results under sym- 
metry assumptions). Hence, property II holds in general for 
the DNMR equations; properties I, III, and I'V hold if shear 
viscosity and heat conduction are absent (with or without 
dynamical coupling to Einstein’s equations); and property 
IT holds with all viscous fluxes present but in the absence of 
chemical potential [75] (with or without dynamical cou- 
pling to Einstein’s equations). Most importantly, the 
DNMR theory has been very successful in phenomeno- 
logical studies of the quark-gluon plasma, particularly in 
numerical simulations of its dynamical behavior; see, e.g., 
Refs. [5,76]. 

We now move to discuss the resumed Baier- 
Romatschke-Son-Starinets-Stephanov theory [77]. In order 
to do so, we need to start with the (plain, not resumed) 
BRSSS theory [77]. This is an effective theory obtained 
from the gradient expansion truncated at second order. As 
such, it involves only the hydrodynamic fields, and the 
equations of motion were chosen in Ref. [77] to be defined 
in the Landau frame. This effective theory-based approach 
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was originally developed for conformal fluids in Ref. [77], 
and the same equations of motion for a conformal system 
were concurrently derived in Ref. [44] through the fluid- 
gravity correspondence, a powerful technique introduced in 
that work which was motivated by the holographic duality 
of string theory [78]. In order to address the issues with 
causality and stability, Baier et al. [77] proposed a MIS-like 
theory with transport coefficients that ensure its agreement 
with the gradient expansion at second order. In the context 
of Ref. [77], this approach provides a resummation of 
higher-order terms and the latter explains the differences 
found, for instance, between rBRSSS and DNMR. 
However, at the linearized level, this resummed BRSSS 
theory shares the same properties of DNMR. Furthermore, 
the techniques used in Ref. [72] can be adapted to establish 
causality for this theory in the nonlinear regime. The local 
well-posedness and hyperbolicity aspects of rBRSSS have 
not yet been established. 

Because the MIS, DNMR, and rBRSSS theories share 
many properties, in particular, the use of extended variables 
that satisfy similar relaxation-type equations, and their 
linearizations about equilibrium agree, they are sometimes 
collectively referred to as Israel-Stewart or Miiller-Israel- 
Stewart theories, Israel-Stewart-like or Miiller-Israel- 
Stewart-like theories, or yet generalized Israel-Stewart or 
Miiller-Israel-Stewart theories. They are sometimes also 
collectively referred to as second-order theories. While 
there is no harm in grouping these theories together in this 
fashion, especially if one is concerned only with their 
general qualitative behavior, it is important to note that 
when it comes to specific features, including properties I- 
IV, the exact form of the equations matters and, therefore, 
the differences among these theories become important. 

The fourth extended theory we would like to briefly 
discuss is the anisotropic hydrodynamics theory [79-83]. 
The latter is, in principle, more general than most 
approaches as it investigates the problem of small devia- 
tions around a given anisotropic nonequilibrium state. 
Formally, this approach involves a resummation in both 
Knudsen and inverse Reynolds numbers, which may be 
interpreted as a generalization of DNMR’s power-counting 
ideas [84]. The equations of motion, which are in practice 
derived using kinetic theory, can be approximated to give 
rise to a MIS-like theory. As such, causality and stability in 
the linearized regime follow from previous results. Nothing 
is known about causality in the nonlinear regime of this 
theory. The local well-posedness and hyperbolicity aspects 
of AHYDRO have not yet been established. 

The above summary highlights how the use of extended 
variables has led to many successes in the study of relativistic 
viscous fluids. These accomplishments seem even more 
impressive when they are contrasted with the fact 
already mentioned that first-order theories (which do not 
employ extended variables) had been largely ruled out for 
decades due to instabilities and lack of causality [25,55]. 
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Such successes nonetheless, it is important to keep in mind 
several actual or potential limitations of the extended theories 
discussed above, as we now discuss. 

First of all, observe that none of the theories MIS, 
DNMR, rBRSSS, or AHYDRO is known to satisfy all the 
properties I-IV. To the extent that they satisfy some of these 
properties, this happens under restrictive assumptions. 
Indeed, in the case of the quark-gluon plasma it is 
abundantly clear that one needs to consider situations 
when all viscous fluxes are present and the chemical 
potential is nonzero [85] (and it is likely that this is also 
true in neutron star mergers [12,14]), in which case none of 
these theories is known to be causal and locally well posed. 
Moreover, while numerical simulations of the dynamics of 
the quark-gluon plasma based on the DNMR equations 
have been carried out for a long time [86-88], only recently, 
with the aforementioned causality results [71,72], one can 
determine regions in the parameter and state spaces for 
which causality holds or fails. When such constraints are 
taken into consideration, it is found that state-of-the-art 
numerical simulations of the quark-gluon plasma violate 
causality [89,90], especially at early times [90]. Although 
further research is required to find out the implications of 
such causality violations to our current understanding of 
those properties of the quark-gluon plasma that have been 
extracted from numerical simulations, such results should 
serve as a definite cautionary tale about running numerical 
simulations of relativistic viscous fluids whose causality 
properties are poorly understood. Furthermore, if causality 
violations can be a real issue in numerical simulations of 
the quark-gluon plasma, which are carried out in flat 
spacetime, the situation is even more precarious in simu- 
lations of general relativistic viscous fluids, such as in 
neutron star mergers. While some simulations have been 
implemented in this setting [11], they rely on a formulation 
for which the key properties I, II, and IV are not known 
to hold. 

Another potential limitation of the extended theories 
discussed above is that they do not seem appropriate for 
describing shock waves [91-93]. This is a potentially 
important limitation given the preponderance of shock 
waves in fluid dynamics, which is aggravated by the recent 
discovery that solutions to MIS-like equations can become 
singular in finite time [94]. Additionally, MIS-like and 
AHYDRO theories are only expected to describe the 
transient regime of dilute gases as their derivation is 
most naturally understood within kinetic theory [36,65]. 
Therefore, their use in other types of systems, such as in 
strongly coupled relativistic fluids, is a priori not justified. 
In fact, it is known that MIS-like equations do not generally 
describe the complex transient regime of holographic 
strongly coupled gauge theories [95-97] (see Ref. [98] 
for the case of higher-derivative corrections). In this aspect, 
we anticipate that the causal and stable first-order theory 
developed here does not describe this transient regime 
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either, despite satisfying properties I-IV. However, this is 
not an issue per se given that the description of such a far- 
from-equilibrium state is certainly beyond the regime of 
applicability of first-order hydrodynamics. 

Finally, MIS-like theories lack the degree of universality 
expected to hold in hydrodynamics as the equations of 
motion themselves change depending on the derivation. For 
instance, the equations of motion in Ref. [77] have different 
terms than in Ref. [65], which is explained by the different 
power-counting scheme employed in those works. This 
situation should be contrasted with theories derived from 
the gradient expansion: although, of course, a plethora of 
different effective theories can be derived in the gradient 
expansion formalism, these different theories can always be 
viewed as particular cases, obtained via different frame 
choices, of the most general expansion truncated at a 
certain order. In fact, an approach of this type is employed 
in this paper; see Sec. IEC. 

Summarizing, despite its undeniable success in advanc- 
ing our understanding of relativistic viscous fluids in 
general, and of the quark-gluon plasma in particular, 
MIS-like and AHYDRO theories still face many chal- 
lenges, especially when it comes to settings where general 
relativity is involved. Thus, it is extremely important to also 
consider alternative theories of relativistic viscous fluids. 
This is especially the case when pursuing the study of 
viscous effects in neutron star mergers [12,14,99,100] and, 
as already mentioned, it is far from clear that the MIS-like 
and AHYDRO approach are the correct approaches for this 
setting. 

In view of the above, it is not surprising that researchers 
have explored other theories of relativistic viscous fluids 
than those discussed so far. A natural place to start such an 
investigation is the gradient expansion, and the simplest 
possibility that includes viscous effects is that of first-order 
theories, i.e., effective fluid descriptions arising as a 
truncation of the gradient expansion at first order. On 
the other hand, since, as said, large classes of first-order 
theories are acausal and unstable, one might naturally 
wonder whether such an approach would be doomed to 
fail. In order to answer this, it is important to understand the 
assumptions involved. While it is true that the acausality 
and instability results [25,64] cover large classes of first- 
order theories, these results apply only to theories that 
satisfy 


U,u,T” = e, (1) 


i.e., only to frame choices that preserve the relation (1). In 
other words, the latter means that an observer moving with 
the fluid always sees the energy density as if it were in 
equilibrium, even for states where entropy is produced. 
Therefore, the construction of stable and causal first-order 
theories remains a distinct possibility as long as one avoids 
constitutive relations that imply (1). First-order theories for 
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which (1) holds are often collectively referred to as the 
(relativistic) Navier-Stokes theory [65], although there is no 
universal agreement on the terminology [35]. 

The physical meaning of (1), as well as of not satisfying 
it, is discussed in Sec. ITC. We also remark that the 
assumptions in Refs. [25,64] imply other special relations 
than (1). But here, for simplicity, we focus only on (1), 
since our goal is not to have a detailed discussion of the 
assumptions involved in those works but rather to illustrate 
how their conclusions apply only for a particular class of 
theories that employ very specific frame choices and, 
therefore, say nothing about first-order theories that employ 
other hydrodynamic frames. In other words, here the reader 
can take (1) as a placeholder for the class of frames that 
are assumed in the instability and acausality results 
[25,64]. Such a class of frames is far from exhaustive. 
Consequently, the results in Refs. [25,64] simply do not 
apply if different constitutive relations are used. 

This motivated researchers to construct stable and causal 
first-order theories of viscous fluids. Important attempts in 
this direction go back to the first decade of this century 
[21,42,47,101]. The first more formal indication that causal 
and stable first-order theories could be constructed if the 
frame choice (1) is avoided is given in Refs. [102-104]. 
These works were also the first ones to carry out a 
systematic study of viscous shocks in relativistic theories, 
a topic that in fact seems to be one of the main goals in 
these references. 

The first construction of a stable and causal first-order 
theory of viscous fluids was carried out by Bemfica et al. 
[32] for the case of conformal fluids (see also Ref. [105] for 
some of the mathematical details of Ref. [32]). These 
results hold with or without dynamical coupling to 
Einstein’s equations. Although Ref. [32] was restricted 
to conformal fluids, it provided an unequivocal proof that 
first-order stable and causal theories are possible, provided 
that one avoids the frame choice (1). Soon thereafter, causal 
and stable first-order theories were obtained by Kovtun [34] 
and by Bemfica et al. [33] for the case of nonconformal 
fluids without a chemical potential [106]—although 
stability was obtained only with the help of a numerical 
investigation, so it might be more precise to say that 
stability was only strongly suggested and not established. 
The resulting first-order theory became known in the 
literature as the BDNK theory [35]. Its local well 
posedness and strong hyperbolicity was established in 
Refs. [107,108]. The stability and causality of the 
BDNK theory in the presence of a chemical potential 
was obtained in Ref. [35] (again, stability in this case was 
inferred only numerically). We also mention the closely 
related results [109,110]. Of course, all these results are 
obtained using frame choices different than (1). Perhaps not 
surprisingly, after these results, the community took a 
renewed interest in first-order theories. See, e.g., 
Refs. [51,98,111—119], and references therein. We remark 
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that choices of frames other than (1) have been studied 
before BDNK in Refs. [42,101,102,120], but, as said, the 
first construction of a stable and causal first-order theory 
was done in Ref. [32] in the case of a conformal fluid. We 
return to the BDNK theory in Sec. II C. In what follows, we 
continue with our brief review of viscous theories. 

Another first-order theory of interest is the Lichnerowicz 
theory [121], introduced in the 1950s but not investigated in 
detail until recently (see references that follow). The 
Lichnerowicz theory has been shown to be causal in the 
(very special) case of irrotational fluids [122] by the second 
author of this paper (see also Ref. [123]). While irrotation- 
ality is too strong of a constraint to be useful for most 
physical applications, Ref. [122] is of interest because it 
initiated the techniques that have since then been employed 
to study the causality of the BDNK theory, including the 
techniques employed in this work. We should also mention 
that the Lichnerowicz theory has found some interesting 
applications in the study of dissipative cosmological 
models [124-127]. 

Another formalism of importance in the study of viscous 
theories is that of divergence-type (DT) theories [128]. In 
this approach, all the conserved quantities describing the 
dynamics of the fluid are obtained from a single generating 
function y which is a function of a dynamical set of 
variables C4 = (€,0,,0,.) (with ¢,, trace-free and sym- 
metric) representing the degrees of freedom of the fluid. For 
example, in the DT approach the energy-momentum tensor 
is obtained as 
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y hed _— P 
aT Og M Og 7) 





DT theories provide a far-reaching subject with many 
important contributions to the physics of fluids, kinetic 
theory, and out-of-equilibrium phenomena. Here, we limit 
ourselves to discuss DT theories with respect to properties 
LIV. See Refs. [2,61,92,128—131] for further discussion of 
DT theories and Refs. [132-134] for applications of DT 
theories to the quark-gluon plasma. 

All information of DT theories is contained in the 
generating function y. Unfortunately, there is no prescription 
on how to construct y, not to speak of how to construct a 
generating function that leads to a theory satisfying proper- 
ties I-IV. In fact, we think it would be more accurate to 
consider the DT approach as a general formalism instead of a 
precisely defined theory or set of theories. That is because 
radically different theories, such as Eckart’s and certain 
types of extended theories, can be cast in divergence-type by 
the choice of a suitable generating function [128]. 

Properties I-IV have been investigated in the context of 
DT theories in Ref. [128]. The authors constructed a DT 
theory that satisfies properties I-IV for states in equilib- 
rium, i.e., when €,=€ Alege Next, they argued that, by 
continuity, these properties will also hold for ¢, sufficiently 
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close to €4|.g- However, no estimate is obtained for how 
close to € Aleg the state €, needs to be. Thus, given any 
nonequilibrium state ¢,, this continuity result does not 
provide any information on whether this specific system 
satisfies the desired properties I-IV. In particular, without a 
quantitative estimate on how small ¢4 — C4|-q needs to be, 
one does not know whether the states €, for which 
properties I-IV hold include states of physical interest. It 
could in principle happen that this continuity argument 
only guarantees the desired properties in a neighborhood of 
Caleg that is orders of magnitude smaller than the size of any 
deviation from equilibrium that one typically considers in 
viscous fluid dynamics. 

Another way of saying this is that the results in Ref. [128] 
are purely qualitative, not providing a quantitative assess- 
ment of their applicability to physical systems. This should 
be contrasted with the precise quantitative results we 
establish here (see Secs. [V—VI) and in the predecessor 
works [32,33,72], which are obtained by employing sub- 
stantially more refined techniques than a general continuity 
argument. In Refs. [92,130-132,134], further results have 
been obtained, but they are all of the same qualitative nature 
as above, relying on precisely the same continuity argument. 
Thus, we believe that a fair assessment of DT theories is that 
they can in principle accommodate properties I-IV, but 
precise conditions ensuring that such properties hold—in 
particular, conditions that allow application to concrete 
physical problems—are yet unknown. 

We finally briefly mention recent formulations of vis- 
cous fluids [135,136] inspired by Carter’s formalism and 
the variational principle [112]. Such formulations address 
some of the properties L-IV but do not establish them in 
completeness. 

Although the review here provided is not exhaustive, we 
believe that it suffices to get across the following main 
point, namely, despite intense work on the subject and 
many different proposals made in the last 80 years, one still 
does not have a theory of relativistic viscous fluids that 
incorporates all relevant viscous fluxes and chemical 
potential while satisfying all the properties I-IV. 
Constructing such a theory is the goal of the present paper. 


C. Summary and discussion of our results 


In this paper, we consider the BDNK theory with 
chemical potential and all relevant viscous fluxes, namely, 
bulk viscosity, shear viscosity, and heat conduction, and 
show that it satisfies all the properties I-IV, i.e., causality, 
stability, local well posedness, and strong hyperbolicity. 
Our results hold in the full nonlinear regime for the fluid 
equations in a fixed background or dynamically coupled to 
Einstein’s equations. We work in 3 + 1 dimensions and do 
not make any symmetry or simplifying assumptions. As 
explained in the previous section, this is the first time that a 
theory of relativistic viscous fluids with all these properties 
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is constructed. In addition, all our results are mathemati- 
cally rigorous and we provide a set of precise inequalities 
among scalar quantities (e.g., shear and bulk viscosity) that 
determine the regions in parameter and state space for 
which properties I-IV hold. Such inequalities are useful for 
numerical simulations as they allow us to check, at each 
time step, whether conditions for causality and stability are 
fulfilled. 

The key conceptual ingredient that allows us to establish 
our results is the realization that the causality and stability 
properties of a theory are intrinsically tied to its hydro- 
dynamic frame. This happens because different choices 
affect the properties of the corresponding partial differential 
equations (PDEs) that describe the evolution of the fluid. In 
particular, we avoid the frame choice (1), which in first- 
order theories leads to acausality and instability. The frame 
choice (1) has a natural intuitive appeal; namely, it states 
that the energy density measured by an observer moving 
with the fluid (1.e., in the fluid’s local rest frame), u,,u,7T"”, 
can be parametrized by a single scalar that can be identified 
with the energy density of the fluid in equilibrium [note that 
(1) holds for an ideal fluid]. It is not surprising, therefore, 
that Eckart and Landau and Lifshitz adopted frames 
satisfying (1). On the other hand, such a simplicity in 
the definition of the hydrodynamic fields out of equilib- 
rium, while desirable, is by no means a fundamental 
property. The key idea underlying the BDNK theory is 
that one should let the fundamental principle of causality 
(and also of stability and local well posedness) dictate 
which frame choices (i.e., parametrizations of T“”) are 
allowed, rather than choose a frame based on nonfunda- 
mental principles and only then investigate properties such 
as causality. In passing, we note that the MIS-like theories 
discussed in this section also adopt (1), although, as just 
said, other frame choices can be made. Different frames 
have been recently investigated in the context of extended 
theories in Refs. [137,138]. 

The idea of exploring different frame choices to con- 
struct a first-order theory that satisfies properties I-IV is not 
entirely new to this work. It was, in fact, the key idea 
employed in the earlier versions of the BDNK theory that 
have been showed to satisfy those properties in some 
particular cases (see Sec. IIB). We next explain what the 
new aspects of this work are, but in order to do so, we need 
to first review some other key ideas employed in the earlier 
constructions of the BDNK theory. 

Since we do not want to make premature frame choices, 
our first step is to consider the most general frame; i.e., we 
write down the most general expression for T”” (and also J¥ 
in the case of the present work since we here consider 
nonzero chemical potential) compatible with the gradient 
expansion truncated at first order; see Eqs. (5) and (6) for 
the precise expression. By considering the most general 
constitutive relations compatible with the symmetries of the 
problem as our starting point, we are in fact applying the 
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basic tenets behind the construction of effective theories 
[48,139-141] to formulate hydrodynamics as a classical 
effective theory that describes the near equilibrium, long 
time, long wavelength behavior of many-body systems in 
terms of the same variables {T,,u’} already present in 
equilibrium. For completeness, we remind the reader that 
an effective theory is constructed to capture the most 
general dynamics among low-energy degrees of freedom 
that is consistent with the assumed symmetries. When this 
procedure is done using an action principle, the action must 
include all possible fields consistent with the underlying 
symmetries up to a given operator dimension and the 
coefficients of this expansion can then be computed from 
the underlying microscopic theory. These coefficients are 
ultimately constrained by general physical principles such 
as unitarity, CPT (charge, parity, and time reversal) 
invariance, and vacuum stability. Analogously, in an 
effective theory formulation of relativistic viscous hydro- 
dynamics, the equations of motion must take into account 
all the possible terms in the constitutive relations up to a 
given order in derivatives that describe deviations from 
equilibrium. The coefficients that appear in this expansion 
can then be computed from the underlying microscopic 
theory (using, for instance, linear response theory [48]), 
being ultimately constrained by general physical principles 
such as causality in the case of relativistic fluids [20] and 
also by the fact that the equilibrium state must be stable; 1.e., 
small disturbances from equilibrium in an interacting (uni- 
tary) many-body system should decrease with time [142]. 
Observe that by considering the most general energy- 
momentum tensor at first order, we are allowing viscous 
corrections to the equilibrium energy density; i.e., one has 


u,uT” = e+ O(T, pn). 


[See Eq. (7) for the precise expression.] Even though this is 
in sharp contrast with (1), in hindsight it seems the natural 
thing to do. After all, it is standard to do precisely the same 
with the pressure, i.e., to split 4A,,T™ into an “equilib- 
rium” part and a “viscous part” (see Sec. IT B) [143]. There 
is no reason not to follow a similar recipe for the energy 
density seen by a comoving observer. 

We next investigate how causality constrains the con- 
stitutive relations. The idea that one should let causality 
determine which frames are allowed in a theory, while 
conceptually powerful, does not tell us how to in practice 
find the appropriate frames. Causality of a theory can be 
determined by computing its characteristics [144]. 
Roughly, the characteristics are hypersurfaces in spacetime 
that correspond to the propagation modes of a theory. For 
example, in the case of Einstein’s equations, the character- 
istics are simply the light cones g,,,v“v” = 0. While in 
principle we can always compute the characteristics of a 
system of PDEs, in practice a brute-force calculation of the 
characteristics seems unattainable for a nonlinear system of 
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PDEs as complex as the BDNK system. In order to be able 
to compute the characteristics, we take a cue from the 
system’s underlying geometric properties. Inspired by 
structures found in the case of ideal fluids by Disconzi 
and Speck in Ref. [145], which need to be recovered in the 
ideal limit, we look for acoustical-metric-like structures. In 
addition, knowing what the characteristics of the system 
should be in some particular limit (e.g., in the conformal 
case that had already been treated) is also helpful to guide 
the calculations. In the case treated here, in particular, we 
already know what needs to be recovered in the limit of 
zero chemical potential. Finally, physical intuition also tells 
us what kinds of modes of propagation should be present in 
the system. In a nutshell, by relying on geometrical and 
physical intuition and an understanding of the causal 
properties of the theory in some particular limits, we can 
have a good educated guess for what the characteristics 
should look like. This allows us to look for a specific 
factorization of the characteristic determinant that points in 
that direction. This is the reason why, in our calculations, 
we group certain terms in certain ways, leading to expres- 
sions that can be managed in the end. Naturally, a brute- 
force approach would not be able to anticipate how one 
should group and factor terms in a way that would allow an 
explicit determination of the characteristics. 

The next step is to carry out a diagonalization of the 
principal part of the equations of motion in order to establish 
strong hyperbolicity. We are able to do so because we have a 
precise understanding of the system’s characteristics. Even 
so, in order to carry out the diagonalization, we need to write 
the system as a system of first-order PDEs (notice that 
V,,T’” = 0 is a system of second-order PDEs because T*” 
involves up to first derivatives of the hydrodynamic fields). 
In doing so, there is the risk of introducing spurious 
characteristics. For example, in the standard linear wave 
equation the characteristics are the light cones. However, 
when one writes it as a first-order system in the standard way, 
the resulting system has a spurious characteristic (it corre- 
sponds, in the language of eigenvalues that can be applied to 
first-order systems, to a zero eigenvalue). While the presence 
of spurious characteristics per se is not an obstacle to 
diagonalization, the more of them there are, the more likely 
there will be obstacles to the diagonalization. Thus, we seek 
to choose as variables for our first-order system quantities 
that have direct physical or geometrical meaning, so that the 
roots of the resulting characteristic polynomial resemble as 
closely as possible the ones of the original system. Of 
course, this does not guarantee diagonalizability. We still 
need to carry out some work mostly technical in nature to 
assure that the system is diagonalizable. But mutilating the 
equations upon rewriting them as first order by introducing 
new, fake features is likely to only make the technical work 
harder or even insurmountable. 

With diagonalization at hand, we can proceed to estab- 
lish local well posedness. The basic idea is that once the 
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system is diagonalized, one can rely on techniques of 
diagonal systems of PDEs. There is a catch, though. The 
diagonalization of the system is at the level of the so-called 
principal symbol (i.e., it is a purely algebraic procedure that 
does not deal directly with differential operators). In order 
to apply it to the actual system of PDEs, one needs to 
introduce pseudo-differential operators, and the quasilinear 
nature of the equations causes further complications as we 
need to deal with pseudo-differential operators with limited 
smoothness. While there are results available in the 
literature for such situations (see, e.g., Ref. [146]), we 
have not found a result that could be directly applied to our 
case. Thus, Bemfica and co-workers developed the neces- 
sary tools in Refs. [107,108] with applications to the 
BDNK equations with zero chemical potential in mind. 
From these techniques and the diagonalization, local well 
posedness follows. 

Finally, let us address stability. For this, one needs to find 
the roots of the polynomial determining the Fourier modes 
of the perturbations. More precisely, only the sign of the 
roots is relevant. Since the corresponding polynomial is of 
high order, there is little hope of determining its roots 
exactly, and even the analysis of the sign of the roots is very 
challenging. Moreover, differently than what happens to 
the causality analysis, geometrical intuition is not of much 
help here because the Fourier modes are not covariant 
quantities. Because of these difficulties, in previous works 
the stability of the BDNK equations was not determined 
rigorously, being obtained numerically or only in the 
homogeneous Lorentz boosted frame [33,35]. Because of 
a new result demonstrated in this paper, this limitation is 
eliminated, as we discuss below. 

We are now ready to discuss specific novelties of the 
present work. While we continue to employ the ideas 
described above and in fact improve on them, especially 
with respect to some of the technical aspects that are more 
challenging for the complete system here considered, we 
want to highlight what are the truly new aspects introduced 
in this work. First, we are able to completely and rigorously 
determine the stability of the system. For this, we rely on a 
new stability theorem, which roughly says that stability in 
the fluid’s local rest frame (which can in general be 
determined because in this case the polynomial for the 
modes simplifies considerably) implies stability in any 
Lorentz boosted frame provided that the system is causal 
and strong hyperbolic; see Sec. VI for the precise assump- 
tions and statement of the theorem. The theorem thus 
establishes a close relationship between causality and 
stability. While connections between causality and stability 
have been discussed before, see Refs. [24,26] and refer- 
ences therein, these results focused on specific theories, 
thus making unclear whether they were due to the specific 
form of the equations of motion or if they were examples of 
a yet undiscovered connection between causality and 
stability as general physical principles. Our theorem, in 
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contrast, is a general theorem that can be applied to many 
different systems, showing that the relationship between 
causality and stability runs deeper and is not a feature of 
specific systems. In fact, we obtain stability of the BDNK 
system by showing that it satisfies the assumptions of the 
general theorem. 

Interestingly, recently, a related theorem was proven in 
Ref. [147], albeit using entirely different methods. The 
results in Ref. [147] also provide further physical intuition 
on the relationship between causality and stability, showing 
that lack of causality allows that dissipation in one Lorentz 
frame be viewed as “antidissipation” (i.e., dissipation 
running “backward in time”) in another Lorentz frame. 
We also note the related work, Ref. [69]. Combined, our 
paper and the works of Refs. [69,147] provide a compre- 
hensive picture of the relationship between causality and 
stability, an idea that was hinted at several times before in 
the literature (see above references) but that had eluded the 
community until now. 

We now discuss strong hyperbolicity. While strong 
hyperbolicity has been obtained for the BDNK theory 
before in the absence of a chemical potential [33,107,108], 
the introduction of a chemical potential causes new severe 
difficulties and the approach used in the case without 
chemical potential does not seem to work. Indeed, in 
Refs. [33,107,108], the choice of variables to write the 
system as first order was based primarily on their physical 
interpretation. For example, the viscous correction to the 
equilibrium energy density was one of the variables chosen. 
As just said, a similar approach does not work here. While 
it is often a good idea to consider variables with a physical 
meaning, the first-order reduction we seek to establish itself 
does not need to carry much physical meaning, so an 
approach employing easily identifiable physical variables 
might not bear any fruit. The first-order system does carry, 
however, some intrinsic geometric properties, such as 
natural decompositions in the directions parallel and 
perpendicular to uv” or the fact that the characteristics of 
the original system are preserved by the reduction to first 
order. Thus, a choice of geometric variables seems more 
appropriate. That is what we have done, considering new 
variables that involve several tensorial decompositions of 
the original variables. This has the extra advantage that 
several tensorial and geometric properties of the fields can 
be used to carry out the difficult calculations needed to 
diagonalize the system. Yet another advantage is that while 
the previous physical choice of variables was specific to the 
form of the BDNK equations, the geometric approach is 
much more general and, thus, can be adapted to other 
theories in that similar tensorial decompositions hold for 
several fluids equations. Therefore, a second novel aspect 
of this work is a new framework to investigate strong 
hyperbolicity in relativistic fluids. We remark that once the 
diagonalization is carried out, we can rely on the techniques 
developed in Refs. [107,108] to establish local well 


021044-12 


FIRST-ORDER GENERAL-RELATIVISTIC VISCOUS FLUID ... 


posedness. Thus, while local well posedness is probably the 
most technical and mathematical aspect of our results, we 
were able to rely more on previous techniques than any 
other of the results we obtain here. 

In addition, it should by no means be overlooked that, 
although the proof of causality provided here follows 
similar ideas as in our earlier work [33], the fact that we 
are now considering the full set of equations makes the 
analysis much more difficult. Thus, a third novelty of our 
work is a substantial improvement of the techniques 
previously employed to analyze causality. From our cau- 
sality analysis, it follows that the characteristics of the 
BDNK theory are the flow lines, sound waves, the so-called 
second sound, corresponding to the propagation of temper- 
ature perturbations [24], and shear waves (plus heat 
diffusion). In addition, when coupling to Einstein’s equa- 
tions is considered, we find another set of characteristics 
corresponding to gravitational waves. 

Finally, as already stressed many times, the main end 
product of this paper is itself a major novelty, namely, the 
first construction of a viscous theory containing all relevant 
fields and satisfying properties I-IV. We accomplish so by 
building and expanding on several previous ideas and also 
by introducing a series of novel ideas, as described above. 

Having discussed the new aspects of our work, we move 
on to discuss how they combine with other aspects of the 
BDNK theory to provide a promising theoretical tool for 
the study of general relativistic viscous phenomena. We 
begin by pointing out that the BDNK theory has been 
shown to be derivable from kinetic theory and holographic 
arguments [32,33,148]. While derivation from kinetic 
theory by itself is not guarantee that a theory is physically 
meaningful since the coarse-grain procedure might intro- 
duce nonphysical features—indeed, recall that the Eckart 
and Landau-Lifshitz theories are derivable from kinetic 
theory—it is reassuring to establish this connection with a 
microscopic theory. As shown in Ref. [148], the derivation 
of BDNK theory from holography can be done in the 
context of the fluid-gravity correspondence [44] by care- 
fully taking into account the presence of zero modes of 
the corresponding differential operators in the holo- 
graphic bulk. 

Next, we should point out that, contrary to MIS-like 
theories, the BDNK theory is capable of handling shocks. 
By this, we mean that Rankine-Hugoniot-type conditions 
can in principle be obtained for the BDNK theory simply 
due to the fact that the BDNK equations are written as the 
conservation laws V,7/” = 0 and V,J“ = 0. Aside from 
this simple observation, viscous shocks have been recently 
studied for the BDNK theory in the case of a conformal 
fluid using numerical methods in Ref. [149], while math- 
ematically rigorously properties were established in 
Ref. [150]. 

At this point, we need to explain the role of shocks in the 
BDNK theory. Since the BDNK theory is an effective 
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theory truncated at first order in the gradient expansion, it is 
expected to be valid when gradients are not very large, 
which is precisely the opposite of shocks. In order to 
explain what we mean by a description of shocks in the 
BDNK formalism, let us consider for a moment an ideal 
fluid. In this case, one also is assuming that gradients are 
small. Alternatively, one may also see this as the limit 
where microscospic length scales are much smaller than the 
length scales associated with the gradients. However, 
shocks are known to develop in solutions of ideal hydro- 
dynamics, and the study of shocks is indeed an important 
topic within the community. To what extent such shocks are 
accurate depictions of the state of the physical system is a 
legitimate question. Nevertheless, once we have decided to 
study shocks in the context of ideal hydrodynamics, the 
formalism allows us to do so in that the equations of motion 
of ideal fluids can accommodate weak solutions (also 
known as distributional solutions) using the Rankine- 
Hugoniot conditions [23]. The same situation happens 
with BDNK: the formalism in principle allows for the 
study of shocks. Whether or not such solutions are 
physical, or accurate in the sense that the results would 
change significantly if the formalism was extended to 
second order, is an important question that is beyond the 
scope of our paper. However, the point we are making is 
that we can, in principle, study shock solutions in the 
BDNK theory. 

In other words, while the derivation of BDNK theory 
rests on the assumption of small gradients, one might try to 
apply it to situations where in principle gradients are not 
small (like shocks), just like it was done before in the 
context of ideal fluids. Although this seems inconsistent, it 
is precisely what it is done when one employs the equations 
of ideal fluids to the study of shocks. Moreover, it is also 
the case that MIS-like theories are often applied to 
situations where gradients are not so small; see, e.g., 
Refs. [84,90,151—155]. It is an intriguing, almost philo- 
sophical, question why one can sometimes still obtain 
meaningful results in such cases, even though shocks are 
formally beyond the regime of validity of any known 
approach to viscous fluids—an important question, how- 
ever, that is beyond our scope here. 

We now discuss another aspect of importance in viscous 
theories, which is entropy production. Naturally, one needs 
the second law of thermodynamics to be satisfied; 1.e., 
entropy production for physically realizable states of the 
system must be non-negative. Before addressing this point 
in the BDNK theory, however, some important points need 
to be highlighted. Strictly speaking, there is no universally 
understood expression for the entropy of a given system out 
of equilibrium, aside from the one given by the Boltzmann 
equation. Thus, while it is useful to define an out-of- 
equilibrium entropy (which must, of course, reduce to the 
definition of equilibrium entropy in the absence of dis- 
sipation), we need to keep in mind that such a definition is 
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not fundamental or even unique. Moreover, the requirement 
that entropy production be non-negative on-shell uncondi- 
tionally, i.e., to all orders in gradients, is certainly too 
stringent. In fact, since a fluid description is an effective 
description, it has a certain limit of applicability. Therefore, 
one should require that entropy production be non-negative 
only within the regime of validity of the theory (which is 
constructed within a certain approximation scheme). This 
point was stressed in Ref. [102] and discussed in detail in 
Ref. [34]. In fact, enforcing non-negative entropy produc- 
tion even in the presence of any size of gradients was part of 
the Eckart and Landau and Lifshitz theories, but the 
resulting theory is unstable and acausal, as seen, showing 
that this requirement by itself is not guaranteed to lead to 
sensible theories in the context of the gradient expansion. 
Non-negative entropy production to all gradients is also a 
guiding principle in the construction of the MIS theory, but 
so far properties I, II, and IV remain open for it. On the 
other hand, the DNMR equations, that are stable, causal (in 
the absence of chemical potential), and are extensively used 
in numerical simulations of the quark-gluon plasma, do not 
have entropy production non-negative to all orders in 
Knudsen and inverse Reynolds numbers, but they should 
have non-negative entropy production within the limit of 
validity of the theory [65]. The same is true for the BDNK 
theory, as pointed out in Ref. [34] and shown in Sec. III A. 
A thorough discussion of the role of entropy in viscous 
theories can be found in Ref. [51]. 

We finally comment on the ability of the BDNK theory 
to describe realistic physical systems. In order to go beyond 
theoretical aspects and make connection with experiments, 
one needs to carry out realistic numerical simulations of the 
BDNK equations. Not surprisingly, given how recent the 
theory is, such investigations are at an initial stage, but the 
results so far have been encouraging. In Ref. [149], the 
authors carry out numerical simulations of the BDNK 
theory in | + 1 dimensions in the case of a conformal fluid 
and compare the results with simulations of MIS ((BRSSS) 
equations in the same setting. They found that for small 
values of the coefficient of shear viscosity, BDNK and MIS 
provide essentially the same evolution, but their dynamics 
differ for larger viscosity values. Given that small viscosity 
is one of the main regimes of interest of both theories 
(higher-order corrections might become relevant in both 
theories if viscosity is not small), this shows that at least in 
this test case the BDNK theory reproduces the well-studied 
and considerably successful behavior of MIS theory. In 
addition, the BDNK theory also reproduces well-known 
behavior considering Bjorken [156] and Gubser [157-159] 
flows, including the presence of a hydrodynamic attractor 
[32]. Further numerical studies of BDNK theory can be 
found in Refs. [160,161]. 

We also stress the obvious point that being a causal, 
stable, and locally well-posed theory are themselves 
fundamental properties that need to be satisfied as a 
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prerequisite for describing actual physical phenomena. 
Thus, while on the one hand a theory possessing these 
properties is only of formal interest if it is not connected to 
experiments, on the other hand, a theory that has some 
phenomenological success but violates, say, causality, 
cannot be taken as an accurate description of real relativ- 
istic physical phenomena. In this regard, we once more 
remark that, in view of the results presented in this paper, 
the BDNK theory is currently the only theory that satisfies 
the fundamental requirements I-III and the additional 
property IV when all viscous contributions and chemical 
potential are incorporated, including in the case when 
dynamical coupling to Einstein’s equations is considered. 


Tl. GENERALIZED NAVIER-STOKES THEORY 


We consider a general-relativistic fluid described by an 
energy-momentum tensor 7“ and a timelike conserved 
current J“ associated with a global U(1) charge that we 
take to represent baryon number. In our approach, the 
equations of relativistic fluid dynamics are given by the 
conservation laws, 

V,J*=0 and V,TY =0, (2) 
which are dynamically coupled to Einstein’s field equa- 
tions: 


R 

Rw = 7 Iw = 82GT,. (3) 
For the sake of completeness, we begin by recalling the 
case of a fluid in local equilibrium [2]. In this limit, one 

uses the following expressions in the conservation laws: 
TY” = ev'u’ + PA” and JF =nu", (4) 
where e is the equilibrium energy density, n is the 
equilibrium baryon density, P = P(e,n) is the thermody- 
namical pressure defined by the equation of state, u“ is a 
normalized timelike vector (i.e., u,u“ = —1) called the 
flow velocity, and A,, = gy, + U,Uy 18 a projector onto the 
space orthogonal to vu“. The thermodynamical quantities in 
equilibrium are connected via the first law of thermody- 
namics ¢-+ P = Ts + yn, where T is the temperature, s is 
the equilibrium entropy density, and y is the chemical 
potential associated with the conserved baryon charge. We 
note that u“V,,e = 0 and u“V,,n = 0 in global equilibrium. 
These are much stronger constraints on the dynamical 
variables than in the case of local equilibrium where, e.g., 
only the combination u“V,e+ (e+ P)V,,u" vanishes. In 
local equilibrium, both u,7“” and J” are proportional to wu” 
and, thus, the flow velocity may be defined using either 

quantity [2]. 

The system of equations (2) and (3) for an ideal fluid 
[defined by Eq. (4)] is causal in the full nonlinear regime. 
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Furthermore, given suitably defined initial data for the 
dynamical variables, solutions for the nonlinear problem 
exist and are unique. The latter properties establish that the 
equations of motion of ideal relativistic fluid dynamics are 
locally well posed in general relativity [16,17]. 

Let us now consider the effects of dissipation. Without 
any loss of generality, one may decompose the current and 
the energy-momentum tensor in terms of an arbitrary 
future-directed unit timelike vector u“ as follows [48]: 


JY = Nu! + J", (5) 


THY = Ey | PAL ul Ov | u’ OF Te (6) 





where N = —u,J*, E€ = u,u,T*, and P = A,,T”/3 are 
Lorentz scalars while the vectors 7” = Aid" S OY 
—u,T’ A», and the traceless symmetric — tensor 
TH = AMAT, with Al, = 5 (AGA, + AGA’ —2AMA,s), 
are all transverse to u,. Observe that this decomposition is 
purely algebraic and simply expresses the fact that a vector 
and a symmetric two-tensor can be decomposed relatively 
to a future-directed unit timelike vector. The physical 
content of the theory is prescribed by relating the several 
components in this decomposition to physical observables, 
which will then evolve [162] according to Eqs. (5) and (6). 

The general decomposition in Eqs. (5) and (6) expresses 
{J",T} in terms of 17 variables {€,N,P,u", J", 
O",T*’\, and the conservation laws in Eq. (2) give five 
equations of motion for these variables. Therefore, addi- 
tional assumptions must be made to properly define the 
evolution of the fluid. As mentioned before, the NS theory, 
including the standard approach in Refs. [15,52], assumes 
that € = e and NV =n. The same assumption is usually 
made in the MIS theory [36], though different prescriptions 
can be easily defined in the context of kinetic theory 
[46,65,163]. A further constraint is usually imposed on the 
transverse vectors, i.e., either 7“ = 0 or QO“ = 0 through- 
out the evolution. For instance, the former gives J“ = nu“ 
and = T” = eutu’ + (P+ TDA” + uO” + u’Qh¥ + T™, 
where II is the bulk viscous pressure (in equilibrium, 
II = 0, QO” = 0, and J’ = 0). In this case, in an extended 
variable approach such as MIS [36], HI, Q”, and JT“ obey 
additional equations of motion that must be specified and 
solved together with the conservation laws, whereas in the 
NS approach these quantities are expressed in terms of uw“, 
€, and its derivatives. 

In this paper, we investigate the problem of viscous 
fluids in general relativity using the BDNK formulation of 
relativistic fluid dynamics. See Secs. I[B and IIC for a 
detailed discussion of the origins of the BDNK theory and 
the conceptual framework that it entails. As explained 
in those sections, the starting point in the formulation of 
the BDNK theory is the most general expression for 
the energy-momentum tensor and the baryon current at 
first order. 
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In practice, the most general expressions for the con- 
stitutive relations that define the quantities in Eqs. (5) and 
(6), truncated to first order in derivatives, are (following the 
notation in Ref. [34]) 





aT 
E=ete- 7 + eVqut" + eyu"Va(u/T), (Ta) 


ueV /T 


P=P+4+n, a + mV u* + 23u°V,(u/T), (7b) 








av 
N=nt VY; : al sr WV qu" oI v3u"V (H/T), (7c) 














AY /T 
Ot = 9, EE ou ,ut + ONY, (u/T), (7A) 
NYY, T 
TH = 1 + Yau Vqul + AMV (u/T), (Te) 
Tw — —2nor”, (7f) 


where of” = AMP V,Ug is the shear tensor. The transport 
parameters {€;,7;,0;,v;,y;} and the shear viscosity 7 are 
functions of T and «. Thermodynamic consistency of the 
equilibrium state (i.e., that ¢, P, and n have the standard 
interpretations of equilibrium quantities connected via 
well-known thermodynamic relations) imposes that y; = 
¥2 and 0, = 0, [34]. The final equations of motion for 
{T,,u*}, which are of second order in derivatives, are 
found by substituting the expressions above in the con- 
servation laws. In the language of Sec. II A, expressions (7) 
for Eqs. (5) and (6) correspond to the most general choice 
of a hydrodynamic frame for a first-order theory. As 
stressed in Ref. [34], it is of course impossible to not 
choose a hydrodynamic frame since the latter actually 
defines the meaning of the variables {7,y,u“} out of 
equilibrium (see Sec. II A for details). 

In fact, in the regime of validity of the first-order theory, 
one may shift {7,v, u“} by adding terms that are of first 
order in derivatives, shifting also the transport parameters 
{€;,7;,9;,V;,7;}, without formally changing the physical 
content of JT” and J“ [34]. However, there are combina- 
tions of the transport parameters that remain invariant under 
these field redefinitions. In fact, the shear viscosity 7 and 
the combination of coefficients that give the bulk viscosity 
€ and charge conductivity o are invariant under first-order 
field redefinitions, as explained in Ref. [34]. Additional 
constraints among the transport parameters appear when 
the underlying theory displays conformal invariance, as 
discussed in detail in Ref. [32] at «= 0, and at finite 
chemical potential in Refs. [34,35] (see also Ref. [110]). 

Hoult and Kovtun [35] investigated Eq. (7) at nonzero 
chemical potential using a class of hydrodynamic frames 
where €3 = 23; = 0; =0. This corresponds to the case 
where there are nonequilibrium corrections to both the 
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conserved current and the heat flux. This choice is useful 
when considering relativistic fluids where the net baryon 
density is not very large, as in high-energy heavy-ion 
collisions. Conditions for causality were derived and limit- 
ing cases were studied that strongly indicated that this 
choice of hydrodynamic frame is stable against small 
disturbances around equilibrium. Further studies are 
needed to better understand the nonlinear features of its 
solutions (well posedness) and also the stability properties 
of this class of hydrodynamic frames at nonzero baryon 
density in a wider class of equilibrium states. 

In this paper, we consider another class of hydrodynamic 
frames that we believe can be more naturally implemented 
in simulations of the baryon-rich matter formed in neutron 
star mergers or in low-energy heavy-ion collisions. Our 
choice for the hydrodynamic frame is closer to Eckart’s as 
we define the flow velocity using the baryon current; i.e., 
J¥ = nu" holds throughout the evolution (y; = v; = 0). 
Clearly, this limits the domain of applicability of the theory 
to problems where there are many more baryons than 
antibaryons so the net baryon charge is large. 

In this case, it is more convenient to use € and n as 
dynamical variables instead of T and «/T because the most 
general expressions for the Lorentz scalar contributions to 
the constitutive relations involve only linear combinations 
of u“V,,€ and V,,u", given that current conservation implies 
that the replacement u/V,n = —nV,u’ is valid. For sim- 
plicity, we choose to parametrize the out-of-equilibrium 
corrections to the scalars as follows [we note that 0; = 05 
and y; = 72, and in practice, 8 out of the 14 parameters in 
Eq. (7) can be set using first-order field redefinitions [34], 
so one is then left with 7, ¢, o, and three other parameters]: 


E=e+7,[u'Vj,e4+ (e+ P)V,u'4], (8a) 
P= P— Vy + aplu'Vje+ (e+ PV yu], (8b) 


where t, and tp have dimensions of a relaxation time and € 
is the bulk viscosity transport coefficient. When evaluated 
on the solutions of the equations of motion, one can see that 
these quantities assume their standard form as in Eckart’s 
theory up to second order in derivatives because € ~ 
e+ O(0*) and P = P-CV,w + O(8") on shell (we 
follow traditional terminology where a given quantity is 
said to be on shell when it is evaluated using the solutions 
to the equations of motion). 

In fact, we remind the reader that in Eckart’s theory [52] 
the energy-momentum tensor is given by ave = 
Eu Uy + (P > CVU") Aw = 2N6 yy + u,Q, + Uy Qy; with 
heat flux Q, = —xT(u’V,u, + AAV,T/T), where x= 
(e+ P)*o/(n’T) is the thermal conductivity coefficient. 
However, as remarked in Ref. [34], in the domain of validity 
of the first-order theory one may rewrite the Eckart expres- 
sion for the heat flux as QO, = oT|(e + P)/nJA‘V,(u/T) 
plus second-order terms. This is done by noticing that 
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(e+ P)wV,u" + AM@V,P =04 O(0") on shell, which 
implies that one may write, using the standard thermody- 
namic relation [(dP)/(e+P)|=[(dT)/T]+[(nT)/(e+P)]x 
d(u/T), 


AYV,T  — nT 
T e+P 





uw V ut t AYV,(u/T) + O(0"). (9) 
Therefore, one can always choose the coefficients such 
that the heat flux OQ” has the same physical content of 
Eckart’s theory plus terms that are of second order on shell. 
We use this to write this quantity as 


O _ gp EHP) 
v n 


+ to[(e + P)u*Vu, + ALV,P], (10) 


AV, (u/ T) 


where Tg has dimensions of a relaxation time. 

In this work, we make the following choice for the 
constitutive relations that give the energy-momentum 
tensor and the baryon current: 














Jt = nu’, (11a) 
T! = (e+ A)utu’ + (P+ TDA” — 270" 
+ ul Q” + uw’ OK, (11b) 
A = 1,[u’V,e+ (e+ P)V,u4], (11¢) 
Il = —CV ut + tplu'Vie+ (e+ P)V,u’], (11d) 
QO” = to(e+ P)u'V,u" + B.AYMV, E+ B,AYMVin, (11e) 
where 
7 OP oT (e+ P) (O(u/T) 
p= rolGe) te (mae) 28 
_ OP oT(e +P) (O(u/T) 
Pn= To (7) n ( on /, aa 


and t,, Tp, and Tg quantify the magnitude of second-order 
corrections to the out-of-equilibrium contributions to the 
energy-momentum tensor given by the energy density 
correction A, the bulk viscous pressure II, and the heat 
flux Q*. In other words, Eqs. (11) and (12) correspond to 
the frame we consider in this work; thus they provide a 
definition of what we mean by the nonequilibrium hydro- 
dynamic fields. 

The reason for considering the constitutive relations (11) 
and (12) is that they lead to a theory satisfying properties I- 
IV, as we show below. We refer the reader to Sec. IT C for a 
discussion of the ideas and techniques that led to the 
particular choice of Eqs. (11) and (12). 
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The equations of motion for the fluid variables are 
obtained from the conservation laws and they can be 
written explicitly as 

WV n+nV,u = 0, (13a) 


wVje+(e+P)V,uw =-u'V,A-(A+T1)V,u4 
—V,, OH — u'V ju, + 2n0,,0%, 


(13b) 
(e+ P)wV uo + AYV{P 
=-(A+T)u’V uv! — APY, I + ALY, (2no") 
4 
— uV,0% - 3 VQ? - 0,0 -O,0, (13c) 


w = 5 (AAV,u, — A2V,u,) is the kinematic 
vorticity tensor [2]. The equations above show that, 
on shell, A ~ 0 + O(07), I~ -C€V,,u" + O(07), and Q, = 
oT [(e + P)/n]AtV,(u/T) + O(0?). Equations (11)(13) 
define a causal and stable generalization of Eckart’s theory 
that is fully compatible with general relativity, as we prove 
in the next sections. We remark that when one neglects the 
effects of a conserved current altogether, the theory reduces 
to the case studied in Refs. [33,34]. For additional 
discussion about the case without a chemical potential, 
including far-from-equilibrium behavior and also the pres- 
ence of analytical solutions, see Refs. [111,116,117]. 


where @ 


A. Entropy production 


It is instructive to investigate how the second law of 
thermodynamics is obeyed in this general first-order 
approach. This was discussed in detail by Kovtun in 
Ref. [34] and, more recently, by other authors in Ref. [51]. 

The standard covariant definition of the entropy current 
based on the first law of thermodynamics TS“ = 
Put —u,T” — pJ*" [36], together with Eq. (11), can be 
used to show that the entropy density measured by a 
comoving observer is given by 


A 
—u,S" =s+ T° (14) 
Note that in our system one finds that A = 0 + O(07) on 
shell. Furthermore, using Eqs. (11) and (13) one finds that 
the divergence of the entropy current is given by 








2no,,0'" TI n 
Vigie eo = _ Ves MV, (u/T 
uS T T pu pale p Q v ale/ ) 
Oo A AAV ,P Auw’V,T 
ea ee | 1 
ii c Meee P| TT m) 


It is crucial to note [34] that in a first-order approach V,,S“ 
can only be correctly determined up to second order in 
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derivatives [recall that in this argument terms such as 
V,,V.¢ and (V,,)(V_@), for any field @, count as second- 
order terms; see Sec. II A]. This means that not all the terms 
in Eq. (15) actually contribute to this expression at second 
order. For instance, when evaluating Eq. (15) on shell one 
must keep in mind that the last two terms in Eq. (15) are 
already at least of third order and must, thus, be dropped. 
A similar argument can be used to show that the term 
TIV,,u" = —¢(V,,w")? + O(03). Therefore, one can see that 





V SH = 2760" ¢(V,uH)? 
B T T 
+ oT (AAV ,(u/T)|[A“°V a(u/T)] + O(0%), (16) 


which is non-negative when n, ¢, o > 0. Hence, there are no 
violations of the second law of thermodynamics in the 
domain of validity of the first-order theory—higher-order 
derivative terms O(0?) in the entropy production can only 
be understood by considering terms of higher order in 
derivatives in the constitutive relations in 7” and J“, which 
is beyond the scope of the first-order approach. 


IV. CAUSALITY 


In order to determine the conditions under which 
causality holds in this theory, we need to understand the 
system’s characteristics. Our system is a mixed first- 
second-order system of PDEs. While the principal part 
and characteristics of systems of this form can be inves- 
tigated using Leray’s theory [21,105,164], here it is simpler 
to transform our equations into a system where all 
equations are of second order. We thus apply u“V,, on 
Eq. (13a). In this case, the conservation laws (2) coupled to 
Einstein’s equations (3) written in harmonic gauge, 
Gf’ Tiy = 9, read 


uP u%7 gn + ndg uO. gu? + B,(n,u,g)0°g = B, (On, Ou, 09), 


(17a) 
(ceu%uP + BAY) gE + B, AV OL gn 
+ p(t. + 19) uy Fu’ + By(e,n, u, g)07g 
= B,(0e, On, Ou, 0g), (17b) 


(Be + tp) UMMA ge + yu AMA? gn + Cy Ogu 


+ Bk(e,n, u, g)0?g = BS (Oe, On, Ou, Og), (17c) 


GP O2gg'” = Bi (Oe, On, Ou, Og), (17d) 


where Oa = 0,0 (using standard partial derivatives), 
p=(e+P), and A(¢By) = (AgAg + ApgB,)/2. The remain- 
ing notation is as follows. We use 0°¢ to indicate that a 
term depends on at most 7 derivatives of @. A term of the 
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form B(0"'), ..., 0% b,)0°¢;, 1 € (1, ..., k}, indicates an 
expression that is linear in 0° #; with coefficients depending 
on at most @, derivatives of d,, ..., 7, derivatives of @;. For 
example, the term (u“0,¢+ 0,u")g%0%,9,5 would be 
written as B(Je, Ou, g)07g (a term of this form is not 
present in our system; we write it here only for illustration). 
The terms B above are top order in derivatives of g and thus 
belong to the principal part, although, as we will see, their 
explicit form is not needed for our argument, whereas the B 
terms are lower order and do not contribute to the principal 
part. We have also defined 


CF = (sop == 2) AMS) + (ptgu*uh —nd®)5f. (18) 


We notice that by taking u“V, of Eq. (13a) we are not 
introducing new characteristics in the system. This can be 
viewed from the characteristic determinant computed 
below which contains an overall factor of wg, to a power 
greater than one. Theorem I below establishes necessary 
and sufficient conditions for causality to hold in our system 
of equations. We show that the assumptions of Theorem I 
are not empty in Sec. VII A. Throughout this paper, we use 
the following definition for the speed of sound c,: 


OP OP n (OP 
2. = | 
4=(%),-(3),+7@) 


where 5 is the equilibrium entropy per particle. Also, we 
define 


=22(2W/D) er WT) pO) | 


. (20) 








Ks 


Theorem I.—Let (€,n,u“, gag) be a solution to Eqs. (3) 
and (13), with w“u,, = —1, defined in a globally hyperbolic 
spacetime (M, gag). Assume that Assumption | 


(Al) p=e+P,t,tg,tp>0 and n,¢,o20. 


Then, causality holds for (¢,n, u", gag) if, and only if, the 
following conditions are satisfied: 








pto > 0, (21a) 
4 2 
[r. (veto +E ; ox, preto 
2 4n 
> 4pt,T9 | tp(pcstg + oK,) — Be ta: >0, (21b) 
2 4n 
2pT.Tg > Te| pestg + F 4 3 tok,)+ptptg>0,  (2Ic) 
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4 
PTT + OKsTp > T, (ncteo ++ a + ox, 
7] 4n 
+ prptg(l—<,) 26 + 3) (21d) 


The same result holds true for Eqs. (13) if the metric is not 
dynamical. 

Proof.—The proof can be reduced to a computation of 
the characteristics of Eq. (17) [164]. Technical details are 
found in Appendix A. 


V. STRONG HYPERBOLICITY 
AND LOCAL WELL POSEDNESS 


In this section, we investigate the initial-value problem 
for Eqs. (3) and (13). The goal is to show that the system is 
causal and locally well posed under very general condi- 
tions. First, we briefly discuss the initial data required to 
solve the system of equations. Then, we rewrite our system 
as a first-order system. We show that this first-order system 
is diagonalizable in the sense of Proposition I. This means, 
in particular, that the system is strong hyperbolic according 
to the usual definition of the term, as in, e.g., Refs. [2,23]. 
The importance of having strongly hyperbolic equations is 
due to its implications for the initial-value problem. As 
already mentioned, one is generally interested in evolution 
equations that are locally well posed [165]. For equations 
with constant coefficients, local well posedness is equiv- 
alent to strong hyperbolicity [166]. For nonconstant coef- 
ficients and nonlinear systems, such an equivalence does 
not hold [167-169]. However, there remains a close 
connection between strong hyperbolicity and local well 
posedness. For most reasonable systems, once diagonaliz- 
ability is available, one can use known techniques to derive 
energy estimates which, in turn, can be used to prove local 
well posedness; see Sec. I1C for more discussion on the 
techniques involved. This is precisely the case for our 
system of equations. Even though our equations consist of a 
system of second-order PDEs, we can use the diagonalized 
system of first-order equations to derive energy estimates. 
Once these estimates are available, we use a standard 
approximation argument as in Refs. [17,170] to obtain local 
well posedness (see Theorem II). 


A. Initial data 


Equations (13) are second order in ¢, n, and wu". 
Thus, initial data along a noncharacteristic hypersurface 
consist of the values of ¢, n, uv“ and their first-order time 
derivatives. Clearly, the initial u“ has to satisfy u“u, = —1. 
Also, it is important to note that Eq. (13a) is first order and, 
thus, the initial data cannot be arbitrary but must satisfy a 
compatibility condition ensuring that Eq. (13a) holds at 
t = 0. Therefore, one can use Eq. (13a) to write the time 
derivative of n in terms of the time derivative of u/ 
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(this feature would also appear in Navier-Stokes theory in 
the Eckart hydrodynamic frame). 

A natural choice to determine the initial conditions for 
the matter sector is to set an initial state that is within the 
regime of validity of the first-order theory and closely 
reproduces Eckart’s theory. First, one can directly extract n 
and uw“ from J” at the initial spacelike hypersurface. Then, 
one sets the nonequilibrium correction to the energy density 
A in Eq. (11) to zero in the initial state, so then the initial 
value for e equals 7T“”u,u, and the first-order time 
derivative of e¢ is defined in terms of the first-order time 
derivative of the flow velocity (plus spatial derivatives that 
are known in the initial state). Clearly, A will be different 
than zero later during the actual evolution, and its value can 
be used to check if the simulations remain within the 
regime of validity of the first-order approach (i.e., |A|/é 
must remain less than unity). Finally, the time derivative of 
the flow velocity can be set by imposing that the second- 
order on-shell term (e+ P)uV,u’ + A“V,P. vanishes. 
Hence, one can obtain the time derivative of the flow 
velocity and all the other required initial data in the regime 
of validity of the first-order approach, emulating Eckart’s 
theory as much as possible. 

We recall that the initial data for the gravitational sector 
has to further satisfy the well-known Einstein constraint 
equations. We briefly make some comments on this in 
Sec. VIII. 





B. Diagonalization and eigenvectors 


In this section, we write Eqs. (3) and (13) as a first-order 
system, as discussed above. For this, we begin defining the 
variables V = u%0,e€, VY = A¥*0,€, W = u%0,n, WH = 
A¥70,n, S# = u°V ,u", Sh = ATV GU", Fy = U7 Oq9,u» and 
vane = Vals By Ae Then, the equations of motion can be 
cast as 


TeU"OgV alr TQ pO, S” le Tepu"OgSy + B.O,V” +P,O0,W =T1; 


(22a) 

tpA“0,V + Tgpu%0,S" + B.u%0gV4 + B,u%0qgW* 
+ nlty’*0,8% = 14, (22b) 
u70,VE — MM@O,V = rf, (22c) 
U"O0qgW# + nAM*0,SY = 14; (22d) 


u20,S% — A%0,S" — X¥440,F , — YAO, Fo = r¥,, (220) 


UOgF 4 = Agr = T6A> (22f) 
MOF] =A OF =F (22g) 
U"Ogé = 1, (22h) 
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u*O Nn = Po, (22i) 

WO. = Figs (22j) 

U"Oq9a = T1As (22k) 

where the r’s are functions of the fields ¢, u”,..., F4, but 


pV 

not its derivatives and A = of for o > f; 1.e., A takes the 10 
independent values 00,01,02,03,11,12,13,22,23,33 with 
repeated index A summing from 00 to 33, 


2 
TH” = —n(AM 62 + AMSf) + (or —o+ 2) AM"6), (23a) 


=op)a 1 
a B) _ 5 [oa” ye — ulead) pa _ ul? AP” AS] (2 ~ 54); 
(23b) 


VA\=Op)a 1 Oo a Sv 
VS B) = 5ul uP) AV55(2 — 54). 


(23c) 
By 6,4 we mean the Kronecker delta in the sense that when 
A = of, then 6, = 6,5, which equals one when o = f# and 
zero otherwise. Also, the terms r may be functions of the 95 
variables. Equations (22) were obtained as follows. 
Equations (22a) and (22b) come from the conservation 
law V,T"” = 0 when projected into the directions parallel 
and perpendicular to u’, respectively. Equations (22c), 
(22d), (22e), and (22g) correspond, respectively, to the 
identities VaV pé = VpVae = 0, VaV pn = VpVan = 0, 
VaV git” —V pV ql’ = Regu? = (Dal Gq —Opl ac)Ue+ terms 
of order zero in derivatives, and 0,09, — O0a9y = 9; 
all contracted with utAe . Equations (22f) is the Einstein 
equation in the harmonic gauge, ie., 0,029, = terms of 
lower order in derivatives, while Eqs. (22h)-(22k) are the 
definitions of V, W (also using the identity u*V jn + 
nV ,u% = W +nS*% = 0 to eliminate W thoroughly), S”, 
and F,4, respectively. We may now define the 95 x 1 
column vectors ¥ and 8 as 


Win 
Y= | w, |. (24) 
Wa 
and B=(r,,...,7y14)’, where yw, =(V,S’,VW’,W, 


Sb, St, S458)" ER”, wo=(F4,F9,F 4, Fa Fa) ER”, 
and wy = (e,n,u’,g,)’ €R'®, to write the quasilinear 
first-order system (22) in matrix form as 


AO, V = B, (25) 


where, here, &* = A® @ ul,6 (@ being the direct sum). 
The matrix A® is split in the following way: 
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A@ —[% 
At = “). (26) 
Osox29 AG 
where 
TU" pte? — BOF BOF pteu"6? preu"5| pt.u%  preu"6s 
tpAlM™ proutdy Puts Buty TO mye yet yea 
=A" 04.4 u* 8, 04.4 04.4 4x4 0434 4x4 
A® = O41 4x4 Oax4 urs) nhKesy ndM6, ndM6; nAM76) (27) 
s\n rr 
Osxr  —AT — Oaxa Osx Osx u"dy O4x4 Osx 
O41 ASE Ox 04.4 0454 0434 udp 044 
O41 AGEL 04.4 04.4 0454 Ver Oana oo, 
| 
while Proof.—The proof of this proposition is very lengthy and 
we refer the interested reader to check all the details and the 
Uy A —-AMZ, AT —AT,o proof presented in Appendix B. 
—A%] u%l 0 0 0 
ae pi 0 " id . ee C. Local well posedness 
g— | Uu : i é P 
¢ i ‘ ae _ ns aia In this section, we establish the local existence and 
~A™Ti9 — O10x10— Otox1o #“Li0 = Ot0x:10 uniqueness of solutions to the nonlinear equations of 
=A" he Ova Oi Opera u"lio motion in Eqs. (3) and (13). 

(28) We begin by noticing that Eq. (13) used the normaliza- 
tion uu, = —| to project the divergence of T,,, and J“ onto 
the directions parallel and orthogonal to uw”. In order to 

and show that the condition uu, = —1 is propagated by the 
flow, it is more convenient to work directly with Eqs. (2) 
Oixto O1n10 O1x10 O1x190 © D110 and (3). In order to complete the system, we differentiate 

ee aeae Me cia 

Gea Usa Uni: “Oi. Daas uM Uy, 1 twice in the uv’ direction: 
O4x10 4x10 04x10 04x10 O4x10 WV p[u2V (ua) — 0. (30) 

O4x10 4x10 O4x10 04x10 O4x10 
Lt= ae ayuka yyisa yyiaa ise (29) We also differentiate VJ" = 0 once, as in Sec. IV: 


Aa wAa wAa wAa wAa 
v y 10 y 11 y 12 y 13 

HAa wAa Aa Aa Aa 
Xv; y 20 y 21 y 22 y 23 


Aa Aa wAa Aa Aa 
Xv; 30 V31 39 V33 


We are now ready to establish that, when written as a 
first-order system as above, the equations of motion are 
strongly hyperbolic. In Sec. VITA, we show that the 
assumptions of Proposition I are not empty. 

Proposition I.—Consider the system (22). Assume that 
(A1) with 7 > 0 holds and that Eq. (21) in Theorem I holds 
in strict form, i.e., with > instead of >. Let & be a timelike 
covector. Then, (i) det(2l°&,) 4 0, and (ii) for any spacelike 
vector ¢, the eigenvalue problem (€, + Aé,)2&%R = 0 has 
only real eigenvalues A and a complete set of right 
eigenvectors R. 


ub .(V_JY) = 0. (31) 
Observe that Eqs. (30) and (31) imply that u“u, = —1 and 
V,,J“ = 0 hold at later times if these hold at the initial time. 
The main result of this section can be found below. 
Theorem Il.—Let (ZX, gap. Rap, €.@n, A, u", i") be an 
initial-data set for the system composed of Einstein’s 
equations (2) and V,J" = 0, where T,, and J” are given 
in Eq. (11). Assume that uy, =-l, n>0O [171], and that 
V,,J“ =0 holds for the initial data. Assume (A 1) with 7 > 0 
and suppose that Eqs. (21) of Theorem I hold in strict form 
and that the transport coefficients are analytic functions of 
their arguments. Finally, assume that a, é,n,u" € HN (x) 
and that &,,,@, A, a € HN-!(Z), N >5, where H™ is the 
Sobolev space. Then, there exists a globally hyperbolic 
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development of the initial data. This globally hyperbolic 
development is unique if taken to be the maximum globally 
hyperbolic development of the initial data. 

Proof.—tThe proof is found in Appendix C. 


VI. NEW THEOREM ABOUT LINEAR STABILITY 


Any ordinary fluid must be stable against small devia- 
tions from the thermodynamic equilibrium state [15]. (We 
only consider systems such that the equilibrium state is 
unique and has a finite correlation length. Therefore, in 
principle, our discussion does not apply to systems where 
the correlation length in equilibrium can become arbitrarily 
large, such as at a critical point.) We recall that in 
equilibrium £, = u,/T must be a Killing vector, ie., 
Vb, + VB, =9, and also V,(u/T) = 0 [36,172,173]. 
In Minkowski spacetime, nonrotating equilibrium corre- 
sponds to a class of states with constant T and mw and 
background flow velocity u“ = y(1,v) defined by a con- 


stant subluminal three-velocity v, where y = 1/V1 —v’. 
(In this paper, we neglect the constant thermal vorticity 
term; see Ref. [172] for a nice discussion of its physical 
content and consequences.) In the local rest frame v = 0 
and the background flow is simply uv“ = (1,0,0,0). In a 
stable theory, small disturbances from the general equilib- 
rium state T > T + 6T(t,x), w> w+ du(t,x), and u“ > 
u“ + du"(t, x) (with u,du“ = 0) lead to small variations in 
the energy-momentum tensor and current, 67"”(t,x) and 
6J"(t,x), which decay with time. 

The standard theories from Eckart and Landau-Lifshitz 
are unstable, as shown by Hiscock and Lindblom many 
years ago [25]. This instability appears because such 
theories possess exponentially growing, hence unstable, 
nonhydrodynamic modes, which spoil linear stability 
around equilibrium even at vanishing wave number. 
(The frequency of a hydrodynamic mode, such as a sound 
wave, vanishes in a spatially uniform state. On the other 
hand, a nonhydrodynamic mode correspond to a collective 
excitation that possesses nonzero frequency even at zero 
wave number.) For Landau-Lifshitz theory at zero chemical 
potential, this instability is only observed when considering 
a general equilibrium state with nonzero v [25,26,73], 
while in the case of Eckart the instability already appears 
even when v = 0. The lack of causality in these approaches 
implies that it is not sufficient to investigate only the static 
v = Ocase in order to determine the stability properties of a 
general equilibrium state where v #4 0, even though such 
states are in principle connected via a simple Lorentz 
transformation. 

The necessity to investigate the stability properties of 
general equilibrium states where v 4 0 makes linear sta- 
bility analyses of viscous hydrodynamic theories very 
complicated. Already in the local rest frame, finding 
whether the linear modes of the system are stable requires 
determining the sign of the imaginary part of the roots of a 
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high-order polynomial, which becomes a daunting task 
when v # 0 (see Refs. [35,70] for recent examples of how 
complicated a v 4 0 analysis can become in BDNK and 
MIS theory, respectively). 

We prove below a new theorem that gives sufficient 
conditions for causal fluid dynamic equations to be linearly 
stable against disturbances of a general nonrotating equi- 
librium state with arbitrary background velocity. In this 
case, proving stability for the local rest frame implies 
stability in any other frame (note that the word frame here is 
used in the standard context of special relativity, i.e., to 
refer to an inertial observer, and has nothing to do with the 
concept of a hydrodynamic frame discussed in previous 
sections, which concerned the definition of hydrodynamic 
variables out of equilibrium) connected to the local rest 
frame via a Lorentz transformation. This general feature is 
expected to hold in any interacting relativistic system; 1.e., 
no issues should appear if one simply observes a given 
system in another inertial frame. We then use this theorem 
in Sec. VII to find conditions under which the hydro- 
dynamic theory presented here is stable. We remark that our 
results can be used to establish stability at nonzero v 4 0 in 
other theories as well, e.g., MIS, as long as the conditions 
discussed below are fulfilled. 


A. Transforming a second-order system of linear 
differential equations into a first-order one 


We begin by showing how one may convert a system of 
linear second-order PDEs into a first order one, as this is 
needed for the theory discussed in this paper. Let the system 
of linearized second-order PDEs be given by 


SMA) g5y>(X) = N(B5y)*, (32) 
b 


where a and brun from | to n, (A)? are differential linear 
operators of order 2, 9(OW) are linear terms containing 
derivatives of the perturbed fields 5 up to order 1, and 
dy'!(X),...,dw"(X) are the perturbed fields (for instance, 
de, On, etc.). We suppose that Eq. (32) arises from the 
conservation laws —u,0,6T°% = 0, AGO,dT% = 0, and 
O,5I* = —ugu?0,6J° + A%A,5J; = 0, where the first 
two come from 0,67” = 0, while the last equation appears 
only when J is included. In this manner, the derivatives in 
the equations of motion in Eq. (32) shall always appear as 
combinations of u%0,, and A*’Q, only. Thus, if the system 
in Eq. (32) has one or more second-order equations, 
it can be rewritten as a first-order system in the N = 5n 
new variables éy“(X) = u70,w"(X) and dyi(X) = 
A‘O,yw“(X). These definitions automatically lead (32) to 
n first-order linear equations. One then needs to supplement 
those with the 4n dynamical equations that are missing. By 
means of the identity O,0g,y“(X) — 0,0,w“(X) = 0, one 
may find the extra 4n dynamical equations 
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u%d, 0974 (X) — AZO, oy"(X) = 0, giving the needed Sn 
first-order dynamical equations, as required. In matrix 
form it becomes 














A%0,5¥(X) + BS¥(X) = 0, (33) 











where A® and are NxN constant real matrices 
and S5Y(X) is a Nx1 column vector with entries 
oy, dy}, ..., dy", 6". This ends the _ procedure. 
However, if one of the equations in Eq. (32) is already 
of first order but contains variables that have second-order 
derivative in other equations, then one can eliminate this 
equation by using it as a constraint to eliminate one of the 
variables. For example, consider the case of the ideal 
current J“ = nu". In this case, the conservation equation 
O,J* =0 becomes u%0,6n(X) + nd,6u2(X) = 0. If TH” 
has shear or bulk contributions, for example, then the other 
equations must have second-order derivatives of du“. Thus, 
one must write 0,6/% = Oas dy + ndyy, = 0, where dy) = 
AZO,u" and dy = u*Ogn. This is a zeroth-order equation in 
the new variables and, therefore, is just a constraint. One 
may use this constraint in order to eliminate the variable 6y 
in the other dynamical equations. Then, in this case one 
ends up with 52 — 1 dynamical equations for the 5n — 1 
fields. 

Finally, we remark that other approaches to viscous 
relativistic fluids, such as MIS, are already written in the 
format (33) in the linearized regime so the procedure to 
reduce the order of the equations of motion described above 
is not needed and one can move directly to the part below. 





B. New linear stability theorem 


To study linear stability, let us expand the perturbed 
fields in the Fourier modes K“ = (iI, k') by substituting 
oY (X) > exp(iK,,X")5¥(K) = exp(I'r + ik;x')6¥(K) in 
Eq. (33). The result is 














iK,A“6¥(K) + BOP(K) = 0. (34) 
Since K* appears, as aforementioned, as combinations of 
—u"K, = (iT — kv!) and A”’K, K, = (w'K,)? +1? +B, 
where k? = k;k', then the direction of k’ is not relevant 
once one keeps v! arbitrary. Thus, we may write 
K* = —n"n,K” + C¥C,K”, where n, is timelike and ¢,, is 
spacelike, with n“n, = —1, n,¢¥ =0, and ¢,¢" = 1 [for 
example, it is common to choose K" = (K°, k,0, 0) so that 
n, and ¢, are (—1,0,0,0) and (0,1,0,0), respectively]. In 
this case we define Q =n,K® and k= €,K% such that 


K¥ = —Qn" + xC¥ [70]. Then, Eq. (34) can be written as 
iQ(—n,A*)6¥(K) = -ix,A%OP(K) — BoY(K). (35) 


The general form of the covectors n and € is ny= 
yn(—1,c') for any c’ such that 0 < cic; < 1 and where 
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Yn = 1/s/1—cle; > 1, and £, = y¢(—d/c;, d') > 1, where 
did; =1 for an arbitrary unitary d' and Y= 


1/\/1 — (dic;)* > 1. From the Cauchy-Schwarz inequality 
(dic;)* < |c'|? (here, |c'| = \/c‘c;), then one obtains that 





Stability demands that the perturbed modes T = I°(k') are 
such that "rp < 0. Now, consider the eigenvalue problem, 


(Ang + C,)A%x = 0, (37) 


where here A is the eigenvalue associated with the right 
eigenvector ¥. 

Proposition IL—If Eq. (33) is causal, then the 
eigenvalues A are real and lie in the range [—1, 1]. 
Furthermore, det(n,A®) # 0. 

Proof.—Causality demands that the roots of Q(é) = 
det(&,A%) = 0 are such that (i) & = €)(&;) € R and that 
(ii) the curves & lie outside or over the light cone. In other 
words, €€, > 0. If one writes €, = An, + €,, where n and 
€ are real, then condition (i) means that A is real. On the 
other hand, since n and € are orthonormal, then condition 
(ii) means that €,6* = —A? + 1 > 0, which demands that 
A? <1, ie., A € [-1,1]. Now, since Q(€) =0 if and 
only if &€ is spacelike or lightlike, this means that 
det(n,A*) # 0. rT 

Theorem IIIl.—Let Eq. (37) have a set of N linearly 
independent (LI) real eigenvectors {r,,..., ty}. If Eq. (33) 
is causal and stable in the local rest frame O, then it is also 
stable in any other Lorentz frame O’ connected to O by a 
Lorentz transformation. 

Proof—The details of the proof are found in 
Appendix D. However, we summarize some steps here. 
Note that causality enables us to invert the matrix (—n,A°). 
Then, it is possible to rewrite Eq. (35) as 


iQ5Y(K)*(R™)-'R-15¥(K) 
= ind5¥(K)' (RT) IR (=n, A°)"'(CA)5U(K) 
— 5P(K)*(RT)“!R“! (=n, A)" BSP(K), (38) 














where the dagger stands for the matrix transpose and 
complex conjugate operations altogether, T stands for matrix 
transpose operation, while R is the square matrix that 
diagonalizes (—n,A*)~!(€,A%), since Eq. (37) has a com- 
plete set of real eigenvectors in R” with only real eigenvalues. 
Then, we can expand 6'¥(K) in terms of these eigenvectors. 
In the proof, it is shown that 6¥(K)'(R7)~'R7'S(K) and 
OP(K)'(R™)-'R7! (—n,A*)7! (€,A%)6¥(K) are real for any 
Lorentz frame. After some work, we demonstrate that, under 
the theorem’s statements, stability reduces to the condition 
that the term 56¥(K)*(R")-!R7!(-n, A®)~'! BOY (K) must be 
greater than or equal to zero. Since this is proven to bea scalar 














021044-22 


FIRST-ORDER GENERAL-RELATIVISTIC VISCOUS FLUID ... 


under Lorentz boosts, it can be computed in any frame. 
Thus, this implies that if the theory is stable in the LRF and 
obeys the other conditions of the theorem, it is stable in any 
other Lorentz frame. 

We note that this result implies that the original system of 
linearized second-order PDEs in Eq. (32) is stable under the 
stated assumptions. 


1, Applying the stability theorem to a toy model 
To illustrate the application of the stability theorem, 


consider the simple model described by the fields @ and y# 
that obey the first-order dynamical linear equations: 


ua — aA," + Ap = 0, (39a) 


uO, — BAFO, = 0. (39b) 


We consider the case where uw is constant [uw = y(1,v'), 
with y = 1/V1—v* and v* = v'v; < 1] as done in the 
stability theorem of the last section. If we write Eq. (39) in 
matrix form as 





A“0,B(X) + 








ix) =0, (40) 





where P(X) = (¢, w’) is a 5 x 1 column vector, 


-| A a 
0451 O44 














and 
ue —ahg 
At = : (41) 
—pAM uw 
are 5 x 5 matrices, the propagation modes w = w(k') are 
obtained by means of the Fourier transform 


W(X) > e!KiX"(K), where K“ =(a,k'), and are the 
roots of det[iK,A* + B] = 0. Let us write w = iT. Then, 
stability requires that Re([) < 0. In the local rest frame, 
these equations are ! = 0 and I? + AV + afk? = 0, where 
k? = k;k'. Then, stability in the LRF implies the conditions 














ap > 0, (42a) 


A= 0, (42b) 


As for the boosted frame obtained by the Lorentz transform 
PT y(P+iv'k;) and ke oI? +h —-y7'(0 + iv‘k;)?, the 
first root is 1 = —iv'k;, which is stable, while the remain- 
ing two roots demand (after a long but straightforward 
computation) 


(43a) 


(43b) 
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To verify stability via the stability theorem proven in this 
paper, we must verify conditions where Eq. (39) is causal 
and if the matrix ®,A% (with ®, = An, + €,, n and € are 
the unitary timelike and spacelike covectors defined in the 
text) has a complete set of eigenvectors in R°. Proposition I 
guarantees that if Eq. (39) is causal, then A € R. In order to 
study causality, we compute the characteristics €, of the 
system, which reduces to the roots of det(A’é,) = 
(u%E,)? [(w'E,)? — aBA"“E,€,| = 0. Causal roots must be 
real and obey €,é" >0, which gives the conditions 
0 < af <1. These conditions, together with stability in 
the LRF, coincide with the conditions obtained by means of 
the above direct calculation. However, if we did not know, 
a priori, the conditions for stability in any frame (which is 
the case when considering higher-order polynomials 
for the modes), we would still have to obtain the eigen- 
vectors of 





u@ a, 
— | -BAHe@, 


a 
—aA v ®, 


,A% 
. ut®,, dh 


(44) 


We can do it firstly by obtaining the eigenvalues A, which 
may be easily obtained by changing €,— ®, in the 
computation of the characteristics. With that result one 
obtains the eigenvalue A that is the root of uD) =0 
with multiplicity 3 and the eigenvalue A) , which give the 
ay — oparn?) po?) = 0. The corre- 
sponding eigenvectors are as follow. 
(i) For ur!) = 0, the system @) aay!) = 0 has as 
eigenvectors the three linearly independent vectors 
given by 





two roots of (uo! 














wef @ 


where {w4}3_, is a set of three linearly independent 


vectors orthogonal to the vector AM7@\, 
(ii) For (wo)? - oparroo = 0, we assume 


ap #0 and obtain the two eigenvectors, 




















a@(2) 
(2) Uu Di 
- | parvo?) | 46) 


Ca 





[Note that in the special case af = 0, the root u*¥, = 0 is 
the only root with multiplicity 5. We end up with two 
distinct situations: first, if a # 0 or 6 4 0 with af = 0, then 
one obtains four LI eigenvectors as can be seen from 
Eqs. (40) and (41). On the other hand, if a = / = 0, then 
the system is already diagonal and the theorem applies 
directly.] Thus, Eq. (46) completes the remaining two 
linearly independent eigenvectors since A, are distinct 
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eigenvalues, giving the five LI eigenvectors. Then, the 
stability theorem states that the system is stable if 2 > 0 and 
0 < af <1 orifA>O0and a= f = 0. Note that there is a 
slight difference from the condition obtained from the 
direct calculation. To wit, it does not include the case af = 
0 with a or # different from zero. The conclusion is that 
stability in any frame does not necessarily imply strong 
hyperbolicity. However, strong hyperbolicity plus causality 
plus stability in the LRF implies stability in any boosted 
frame. In other words, stability may occur outside the 
conditions imposed by the theorem. 


2. Applying the stability theorem to the MIS system 


As another example of the usefulness of Theorem II, 
let us briefly comment how it can be used to recover 
the stability conditions of the MIS equations [24] in the 
presence of bulk viscosity. More precisely, we take the 
MIS-like equations studied in Ref. [71] where only bulk 
viscous effects have been considered. In that case, it was 
proven that there exist conditions such that the system of 
PDEs is nonlinearly causal and symmetric hyperbolic; 
hence the principal part of the equations is diagonalizable. 
The linear version of such equations forms a system that is 
also symmetric hyperbolic and the conditions for stability 
needed for the application of Theorem III can be shown to 
agree with those found in Ref. [24] for the case where only 
bulk viscosity is present. 


VII. CONDITIONS FOR LINEAR STABILITY 


We now apply the theorem proved in the last section 
to determine conditions that ensure the stability of the 
hydrodynamic theory proposed in this paper. Let us first 
define 


D=pcr(t, +79) +64 





+ OK, (47) 








= Hl 35), an),~ (an), ae 


(48) 


where K,=(Tp2/n)[A(u/T)/Ie]s=Ke-+ky, Ke = (Tp*/n)x 
[A(u/T)/Ae],. Ky = (Tp)[O(u/T)/An],, and p.=(OP/2e),. 
Standard thermodynamic identities imply that p/.«, — 
cx, >0, then E>0 from (Al). By assuming the 
Cowling approximation [174] with g, =.= 
diag(—1,1,1,1) and 6g,, =0, we find that the system 
described by Eq. (13) is linearly stable if it is causal within 
the strict form of the inequalities in Eq. (21) together with 
the additional restriction 7 > 0 in (A1) and 
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(t.+tg)|B| = t-tgD = pcrt.t9(t. +79), (49a) 
(t, + Tg) |B|D + pt,to(t, + t9)E 
> t.t9D* + p(t, + tg)°C, (49b) 
c3D —E > pce}(t. + t9), (49c) 
(t. + tg) [|B|(c2D — 2E) + 2c2pt,t9E + CD] 
> 2c3p(t, + t9)°C + t.tgD(c2D — E), (49d) 
|B|D[C(z, + tg) + Et,tg] + 2pt,t9(t, + tg)CE 
> pC*(t, + tg)” + tetg(CD* + pt.tgE’) 
+ B’E(t, +79), (49e) 
where B and C are given by 
B= 2 zl 50 
=-T,| pestg + + 3 +ok;}—ptptg, (50a) 
2 4n 
C= tp(pcstg + oK;) — Be| ¢ + 3) (50b) 


as in Eq. (A5), with |B] = —B > 0 from Eq. (21c) in the 
strict form. 

To prove the statement above, as before we may expand the 
perturbations 6 = (de, du“, dn) in Fourier modes by means 
of the substitution 65¥(X) > exp[T7(It 4+ k;x’)|6P(K), 
where K* = (iI’, k') is dimensionless due to the introduction 
of background temperature T in the exponent. We begin by 
proving stability in the local rest frame, where the modes are 
the roots of the shear and sound polynomials, 


shear channel: ZI? + 7k? +T = 0, (51a) 





sound channel: agl° + a,I++ a,03 +430? +a4l +a5 =0, 





(51b) 
where k* = k'k; and 
a) = 7,79, (52a) 
a; =7, +79, (52b) 
dy =1+k?|Bl, (52c) 
a, = RD, (52d) 
ag = 2k? + KC, (52e) 
as = KE. (52f) 
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We defined the dimensionless quantities tg = T79, 
7, =Tt., 7 =Tn/p, B=(T?/p)B, C=(T?/p)C, D= 
(T/p)D, and E = (T/p)E. From the second inequality in 
Eq. (21c) in its strict form one obtains that B < 0 (see the 
definition of a2). The analysis of stability in the LRF goes as 
follows. 

Shear stability conditions.—The second-order polyno- 
mial (51a) has two roots with 'p < 0 only if tg > 0 and 
n => 0, which is in accordance with assumption (A1). One 
can see that tg clearly acts as a relaxation time (the same 
role is played by the shear relaxation time coefficient 7, 
present in MIS theory) for the shear channel, which 
ensures causality. In fact, the condition tg > 0 is clear 
since the leading contribution to the nonhydrodynamic 
frequency in this channel goes as 1/tg at zero wave 
number. 

Sound stability conditions.—As for the sound channel in 
the rest frame, by means of the Routh-Hurwitz criterion 
[175], the necessary and sufficient conditions for Ip < 0 
are (i) dg, a; > 0, (ii) ajay — aga3 > 0, (iii) a3(ayaz—- 
dy43) — A, (a,a4—dyas) > 0, (iv) (4;a4—agas)[a3(a,az— 
Ay43) — 41 (4144 —aods)] —as(a,ay—aga3)” > 0, and 
(v) as > 0. Condition (i) is already satisfied from (A1). 
Condition (ii) corresponds to the first inequality 
in Eq. (49a), while (iti) is the second inequality in 
Eqs. (49a) and (49b). Condition (iv) corresponds to 
Eqs. (49c)-(49e). Given that E > 0, thus, when E = 0 
and (i)-(iv) are observed, then Ip <0, which is in 
accordance with stability. Also, if k =0, then Tp <0 
(three zero roots and two negative roots) because 
dp, 4,,a, > 0 from (Al). Hence, the system is linearly 
stable in the local rest frame. 

We remark that our system displays three types of 
hydrodynamic modes and three nonhydrodynamic modes. 
In the small k expansion that typically defines the linearized 
hydrodynamic regime, our shear channel gives a diffusive 
hydrodynamic mode with (real) frequency w(k) = 
—ik’n/(e + P) +--+, while in the sound channel one finds 
proper sound waves with w(k) = -+c,k — ik’T,/2+--- 
and also a heat diffusion mode with w(k) = —iDk? +---, 
where D ~ o, andl’, =T,(7,, 0) just as in Eckart theory 
(see Ref. [35] for their detailed expressions). Therefore, our 
theory has the same physical content of Eckart’s theory in 
the hydrodynamic regime. On the other hand, the shear 
channel has a nonhydrodynamic mode with frequency 
given by w(k) = —i/tg +---, while the sound channel 
has two nonhydrodynamic modes with frequency w(k) = 
—i/t, +--+ and w(k) = —i/tg +--- in the low k limit. 
These nonhydrodynamic modes parametrize the UV behav- 
ior of the system in a way that ensures causality and 
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stability, making sure that the theory is well defined 
(though, of course, not accurate) even outside the typical 
domain of validity of hydrodynamics. 

The complete proof of linear stability demands an 
analysis of the linearized system around an equilibrium 
state at nonzero velocity. In this regard, we shall use the 
results presented in Sec. VIB. We first write the system in 
Eq. (13) as a first-order linear system of PDEs. Then, since 
we already have proven causality and also linear stability in 
the LRF, it remains to be shown that the first-order 
counterpart of Eq. (13) is diagonalizable in the sense of 
Eq. (D2). This is done below. 

First-order system.—Following Sec. VIA, we may 
define 6V = u°0,6e, dV" = A¥*0,6e, SW = u%0,6n, 
bYVE = A¥*90,6W, dS" = u%0,0u", SS) = AF0,6u’". 
Since the current is ideal, i.e., J“ = nu”, then the linearized 
conservation equation 0,,dJ“ = 6W + ndS7 = 0 enables us 
to eliminate 6W from the new system of equations. Hence, 
the first-order equations become 


TU" qSV + ptQ,5S* + BeOg5V" + By Oqg5W* 


+ pt,u*%0,6S) + dV + pdS? = 0, (53a) 


TpAH*O,6V a ptgu?d,6S" te B,u*0 dV" + B,U%OgOVVE 


+ 1/7, 6S% + p6V! + pl dW" + p6S“=0, — (53b) 
u70,5V" — AF20,6V = 0, (53c) 

u70,0VW" + nA“*0,6S% = 0, (53d) 

u",6S% — A%0,6S" = 0, (53e) 


where D', = (OP/On), and 
2. 
THs = —n( AM 68 + AMS) + (v1 ee ee =) Aes). (54) 


The supplemental equations (53c)-(53e) come from the 
identities O,0,6€ — 0,0,6€ = 0, O,0,6n — 0,0,6n = 0, 
and 0,0,;6u" — 0,0,6u" = 0, respectively, when con- 
tracted with u*A/. In particular, in Eq. (53d) we have 
substituted 6bW = —ndS/ that comes from the conservation 
equation of J“. Then, we may write Eq. (53) in matrix form, 
A20,08(X) + BY(X) =0, were 6¥(X) is the 29x 1 
column matrix with entries 6V,6S’,dV", d5W’,dS6, 
6S, 5S5, 6S%, 
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T,U% PTQOy B69 B00 

TpAM* = ptgu*dy Bud B,u%dy 

—AKS 044 uaa 

At = O44 O44 O44 urd, 

Osx1  —ABH — Oaxa O44 

O4x1 Ato 04.4 O44 

Osx1  —ASH — Oaxa O44 

O4x1 ASSL Ox 4x4 

and 

L O01x4 O1x4 1x4 poy po, po, po; 
Osx1 po peor Prov O4x4 04x4 O4x4 O44 
O4x1 O4xa O4x4 O44 O4xg Ogxa Odxa Oana 
| O4x1 Oana O4xa O4x4 O4x4 Osxa O4xg Oana 
© [Oar Osa Oag Oded Oded Oad Oded Osa 
O4x1 O4x4 O4x4 O44 O4xg O4x4d Odxa Oana 
O4x1 O4xa O4x4 Ona O4xg Ogxa Odxa Oana 
O4x1 O4xa O4xg Oaxa O4x4 Osx Oana Osx 
(56) 


We must now obtain the eigenvectors of Eq. (37). However, 
note that A® above is exactly the same as the matrix A%, in 
Eq. (27) with the difference that now the coefficients of A® 
are constants. We have already proven in Sec. V that the 
matrix A’, in Eq. (37) has real eigenvalues and a complete 
set of eigenvectors in R*’. The same solution is true for A“ 
in Eq. (37) if we change €, > n, (and also A% > A®) in 
the results for the matter sector in Sec. V. Thus, the 29 x 29 
matrix (—n,A°)¢ rN is diagonalizable, completing the 
requirements from Theorem III. This shows that the theory 
is linearly stable in any other reference frame O’ connected 
via a Lorentz transformation. Therefore, one then obtains 
that our set of linearized second-order PDEs is stable in any 
equilibrium state. 


A. Fulfilling the causality, local well posedness, 
and linear stability conditions 


We now give a simple example that illustrates that the set 
of linear stability conditions (and consequently, causality 
and local well posedness, since those are part of the linear 
stability conditions) is not empty. Let us analyze the case 
where tg = T, and tp = c?t,, assuming an equation of 
state P = P(e), with c? = p’, = 1/2. Also, assume that 
€ + 4n/3 > 0 (their specific values are not relevant as far as 
they are positive and 7 > 0 for the sake of the stability and 
well-posedness theorems). Then, one may easily verify that 
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pt,u°S? pt,u%6! pt,u%62 pt,u%53 
nm Oa 11 aaa 32 
04.4 0, x4 O44 044 
nd¥?59 nAv*6)  nAH*6 — nA#*653 (55) 
u*dy O44 0434 Ope | 
O44 u*dy O4x4 O4x4 
O44 O44 u*y O44 
O4x4 O44 O44 urd, 








the causality conditions (21) hold in their strict form, as 
required, and that the remaining conditions (49) are also 
observed when pt, = 8(¢ + 4n/3), kK. = K,/2 = 1/4, and 
in the three different situations, namely, o/(¢ + 47/3) = 0, 
1/4, and 1. 


VIII. CONCLUSIONS AND OUTLOOK 


In this work, we presented the first generalization of 
relativistic Navier-Stokes theory that simultaneously sat- 
isfies the following properties: the system, with or without 
coupling to Einstein’s equations, is causal, strongly hyper- 
bolic, and locally well posed (see the content of Theorems I 
and II); equilibrium states in flat spacetime are stable 
(consequence of Theorem III); all dissipative contributions 
(shear viscosity, bulk viscosity, and heat conductivity) are 
included; and finally the effects from nonzero baryon 
number are also taken into account. All of the above hold 
without any simplifying symmetry assumptions and are 
mathematically rigorously established. In addition, entropy 
production is non-negative in the regime of validity of this 
effective theory. 

This is accomplished in a natural way using a general- 
ized Navier-Stokes theory containing only the original 
hydrodynamic variables, which is different than other 
approaches where the space of variables is extended (such 
as in Miiller-Israel-Stewart theory). However, it is impor- 
tant to remark that the meaning of the hydrodynamic 
variables in our work is different than in standard 
approaches, such as Refs. [15,52]. In fact, in the context 
of the formalism put forward by Bemfica et al. [32,33] and 
Kovtun [34], our formulation uses a definition for the 
hydrodynamic variables (i.e., our choice of hydrodynamic 
frame) that is not standard as there are nonzero out-of- 
equilibrium corrections to the energy density and there is 
energy and heat diffusion even at zero baryon density. 
Despite these necessary differences (imposed by causality 
and stability), the theory still provides the simplest causal 
and strongly hyperbolic generalization of Eckart’s original 
theory [52], sharing the same physical properties in the 
hydrodynamic regime (for instance, both theories have 
the same spectrum of hydrodynamic modes). However, 


021044-26 


FIRST-ORDER GENERAL-RELATIVISTIC VISCOUS FLUID ... 


differently than Eckart’s approach, our formulation is fully 
compatible with the postulates of general relativity, and its 
physical content in dynamical settings can be readily 
investigated using numerical relativity simulations. In fact, 
we hope that the framework presented here will provide the 
starting point for future systematic studies of viscous 
phenomena in the presence of strong gravitational fields, 
such as in neutron star mergers. 

Motivated by the task of establishing stability of general 
equilibrium states in flat spacetime, in this work we also 
proved a new general result (see Theorem III) concerning 
the stability of relativistic fluids. In fact, we found con- 
ditions that causal relativistic fluids should satisfy such that 
stability around the static equilibrium state directly implies 
stability in any other equilibrium state at nonzero back- 
ground velocity. Theorem III is very general and its regime 
of applicability goes beyond BDNK theories and it could 
also be relevant when investigating the stability properties 
of other sets of linear equations of motion as well. In this 
regard, see the discussion in Sec. II B, and see also Secs VI 
B 1 and VIB 2 for further examples of the applicability of 
Theorem III. 

Our generalized Navier-Stokes theory can be used to 
understand how matter in general relativity starts to deviate 
from equilibrium. An immediate application is in the 
modeling of viscous effects in neutron star mergers. Our 
approach can be useful in simulations that aim at determin- 
ing the fate of the hypermassive remnant formed after the 
merger of neutron stars, hopefully leading to a better 
quantitative understanding of their evolution and eventual 
gravitational collapse toward a black hole. Differently than 
any other approach in the literature, the new features 
displayed by our formulation and its strongly hyperbolic 
character make it a suitable candidate to be used in such 
simulations. This will be especially relevant also when 
considering how viscous effects may modify the gravita- 
tional wave signals emitted soon after the merger [12,14]. 
In this regard, we remark that previous simulations per- 
formed in Ref. [11] employed a formulation of relativistic 
viscous hydrodynamics where the key properties studied 
here (causality, strong hyperbolicity, and local well posed- 
ness) are not known to hold in the nonlinear regime. 

Our work is applicable in the case of baryon-rich matter, 
such as that formed in neutron star mergers or in low- 
energy heavy-ion collisions. The latter include the exper- 
imental efforts in the beam energy scan program at RHIC 
[176], the STAR fixed-target program [176], the HADES 
experiment at GSI [177], the future FAIR facility at GSI 
[178], and also NICA [179]. For a discussion of viscous 
effects in low-energy heavy-ion collisions at nonzero 
density, see Refs. [85,113,180]. High-energy heavy-ion 
collisions, such as those studied at the LHC, involve a 
different regime than the one considered here where the net 
baryon number can be very small and, thus, that case is 
better understood using a different formulation such as the 
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one proposed in Ref. [35], also in the context of the BDNK 
formalism. 

In our approach, we only take into account first-order 
derivative corrections to the dynamics. Therefore, the 
domain of validity of our theory is currently limited by 
the size of such deviations. Hence, further work is needed 
to extend our analysis, incorporating higher-order deriva- 
tive corrections, to get a better understanding of what 
happens as the system gets farther and farther from 
equilibrium. In this context, it would be interesting to 
extend our equations to include second-order corrections 
and consider also, more generally, the large order behavior 
of the gradient expansion in an arbitrary hydrodynamic 
frame. The latter will be different than most approaches to 
the gradient expansion since in BDNK the constitutive 
relations contain time derivatives even in the local rest 
frame of the fluid. This essential difference has important 
consequences in a kinetic theory formulation; see the 
original references [32,33]. The large order behavior of 
the relativistic gradient series has been recently the focus of 
several works [84,181—194], and it would be interesting to 
extend such analyses to include the type of theories 
investigated here. 

There are a number of ways in which our work could be 
extended or improved. First, it would be useful to obtain a 
better qualitative understanding of why some hydrody- 
namic frames (such as the Landau-Lifshitz frame or the 
Eckart frame) are not compatible with causality and 
stability in the BDNK approach, given that the situation 
is different in other formulations. In fact, the Landau frame 
seems to display no significant issues in the case of MIS- 
like theories in the nonlinear regime at least at zero 
chemical potential, as demonstrated in Ref. [72]. Perhaps 
a more in-depth investigation of how BDNK emerges in 
kinetic theory, going beyond the original work done in 
Refs. [32,33], can be useful in this regard (see also the 
recent work [148]). Also, it would be interesting to use the 
BDNK approach to investigate causality and stability in 
more exotic cases, such as in relativistic superfluids. 
Furthermore, the inclusion of electromagnetic field effects 
in the dynamics of relativistic viscous fluids can also be of 
particular relevance, especially in the context of neutron 
star mergers [195] and high-energy heavy-ion collisions 
[196]. This problem has been recently investigated using 
other formulations of viscous fluid dynamics, see for 
instance Refs. [197-202], and also most recently in the 
BDNK approach in Ref. [203]. Consistent modeling of 
relativistic viscous fluid dynamics coupled to electromag- 
netic fields can also be relevant to determine the importance 
of dissipative processes in the dynamics and radiative 
properties of slowly accreting black holes, as discussed 
in Ref. [197]. 

Further work needs to be done to understand the global 
in-time features of solutions of relativistic viscous fluid 
dynamics. For instance, one may investigate the presence 
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of shocks, which is a topic widely investigated in the 
context of ideal fluids [21,145,204—206] and was done in 
Ref. [94] for the MIS theory (see Sec. ITB for further 
discussion on shocks). The importance of hydrodynamic 
shocks has been recognized both in an astrophysical setting 
[197] as well as in the study of jets in the quark-gluon 
plasma [207-219]. We also remark that one task that we 
have not done here was the construction of initial data for 
the full Einstein plus fluid system by solving the Einstein 
constraint equations. We believe that standard arguments to 
handle the constraints [21] will be applicable in our case. 
This will be investigated in detail in a future work. 

We believe our work will also be relevant to give insight 
into the physics of turbulent fluids embedded in general 
relativity. The fact that the equations of motion of the 
viscous fluid must be hyperbolic in relativity stands in 
sharp contrast to the parabolic nature of the nonrelativistic 
Navier-Stokes equations, usually employed in studies of 
turbulence. Recent works in Refs. [18,220] tackled the 
problem of turbulence in the relativistic regime and our 
formulation may be very useful in this regard, as it provides 
a simple strongly hyperbolic generalization of Eckart’s 
theory that is fully compatible with general relativity. 

In summary, in this paper we propose a new solution to the 
question initiated by Eckart in 1940 concerning the motion of 
viscous fluids in relativity. Our approach is rooted in well- 
known physical principles and solid mathematics, displays a 
number of desired properties, and extends the state of the art 
of the field in a number of ways. Potential applications of the 
formalism presented here spread across a numbers of areas, 
including astrophysics, nuclear physics, cosmology, and 
mathematical physics. This work establishes for the first 
time a common unifying framework, from heavy-ion colli- 
sions to neutron stars, that can be used to discover the novel 
properties displayed by ultradense baryonic matter as it 
evolves in spacetime. 
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APPENDIX A: PROOF OF THEOREM I 


We only consider the 10 independent components of the 
metric and, thus, this system of equations can be written in 
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terms of a 16 x 1 column vector '¥ = (e,n, u’, g,,), and its 
equation of motion in Eq. (17) can be expressed in matrix 
form as IN(O)'¥ = N, where MN contains the G terms that do 
not enter in the principal part. The matrix (0) is given by 


M(O) (2) 
06x10 L109 #O25 J” 


MO) (Al) 


where the 6 x 10 matrix b(Q) contains the B terms and 


0 utuP nd" wy) 
M(Q)= | (cut? +f.A%)  B,A% p(t, +t) ul*op? 
(B.+tp)u AP” Bult AAe cub 
x Oop: (A2) 


The system’s characteristics are obtained by replacing 
Og > €, and determining the roots of det[9%(é)] = 0. 
The system is causal when the solutions for €,= 
(E9(&;), €;) are such that condition 1 (Cond-1) &, is real 
and condition 2 (Cond-2) €,¢“>0 [21]. It is easy to see that 
det[IMN(E)] = (€,E7)!° det[M(E)]. The roots associated with 
the vanishing of the overall factor (€,é%)!° = 0 coming 
from the gravitational sector are clearly causal. The 
remaining roots come from det{M(é)] = 0, which we will 
investigate next. 

We first define b=u%E, and v%= AP Es, which 
gives £,=—bugt+vq and &,€*=—b?+0-v, where 
vev= AE Ey. We proceed by also defining the tensor 


De Od 


= (cop -f- =) vHE, + [ptob? —n(v-v)|&, — (A3) 


which gives 
0 b* nbé, 
det[A(é)] =det | 7.b°+B,(v-v) Bny(v-v) p(te+t9)bv, 
(Be Tv Tp) bot B, bu" a 
=—b*Iprob®—n(o-)} 
x [Ab* + Bb?(v-v) +C(v-v)] (A4a) 
= =preotlwre,) II [Cae ~ Cah PEE s)"*, 


a=1,+ 








(A4b) 
where, to shorten notation in Eq. (A4a) we defined 


A= ptto: (A5a) 


4 
=-t, (pc? 4G =p = + ox; ) — Ptptg, (A5b) 


coniniteren)=a(c%), (so 
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and used the fact that £, + nf,/p = tgc? + ox,/p. In 
Eq. (A4a) it becomes evident that assumption (Al) guar- 
antees that v“ # 0, eliminating one of the possible acausal 
roots. From Eqs. (A4a)-(A4b) we defined n, = 3, n+ = 1, 
c, = [n/(pt,)], and cx = [(—B + VB? — 4AC)/2A]. Note 
that since £*£, = —b* + (v- v), the roots in Eq. (A4b) can 
be cast as b? = c,(v-v). Then, (Cond-1) demands that 
Cg € R together with c, > 0 and (Cond-2) that c, < 1 for 
causality [221], which comes from the fact that the root 
b?>=c,v-v must obey €f/=-b?+v-0= 
(1 -—c,)v-v > 0. Thus, causality is ensured if 0 < c, < 
1 in the matter sector. Clearly, the root b = u%€, = 0 is 
causal. Also, the six roots related to c, are causal when 
Eq. (21a) is observed. As for the roots c.., they are real if 
B? — 4AC > 0, i.e., if the first inequality in Eq. (21b) holds. 
On the other hand, cx > 0 is obtained whenever c_ > 0, 
which is guaranteed if —B >O [second inequality in 
condition (21c)] together with C > 0 [second inequality 
of Eq. (21a)], while c, < 1 is ensured if c, < 1, which 
demands that 2A + B > 0 [first inequality in condition 
(21c)] and A+ B+C > 0 [condition (21d)]. a 

We observe that, although we employed the harmonic 
gauge to calculate the system’s characteristics, the causality 
established in Theorem I does not depend on any gauge 
choices. This follows from well-known properties of 
Einstein’s equations [22] and the geometric invariance of 
the characteristics [144]. See the end of Sec. V C for further 
comments in this direction. 
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The analysis above and the conditions we obtained for 
causality are valid in the full nonlinear regime of the theory. 
However, we remark in passing that the principal part 
concerning only the fluid equations would have exactly 
the same structure if one were to linearize the fluid dynamic 
equations about equilibrium with nonzero flow in 
Minkowski spacetime. This is a generic feature of the 
BDNK approach (at least, when truncated at first order), 
i.e., the analysis of the system’s characteristics, and thus of its 
causality properties, is formally the same in the nonlinear 
regime and in the linearization about a generic equilibrium 
state. This is not, however, a general feature of hydrodynamic 
models as it does not hold in MIS-like theories. In fact, as 
discussed at length in Refs. [71,72], in MIS the thermody- 
namic fluxes explicitly enter in the calculation of the 
characteristics, but they are not present in the linear analysis. 


APPENDIX B: PROOF OF PROPOSITION I 


To prove (i) we may compute the determinant 
det(E,2%) = det(E,A%,) det(E,A%)(u7é,)'®. Note — that 
u°E, #0 if € is timelike. We must then look into the 
matter and gravity sector in what follows. We again define 
b=u%§, and v4 = AM@E,, v-v = AME, ¢,, and introduce 


mH 
BH 


Wha 
= vl Ca 


2 
—n(v- 0) 8 — nvté, + (v= -¢+ 21) vv, (B1) 














to obtain 
| 
Tb ptot, Piece Baby pt-bee pt,b6, pt.b6; pt-b6> 
tev! ptobdt Pb fybo TUE, ME, TE, TS, 
=O Daya bY aya 04,4 O44 04.4 O44 
mu aN) ws sd Hw S3 
det(E,A%) = det Oayy Onyx Oaya bey nud, nvotd, nod, nod, 
Osx1 = 05 Ota Oya be, O4x4 04.4 04x4 
Osx1 018 Oa Oye O44 bey O44 O44 
Osx. v2) axa Ox O4x4 O4x4 be Osx 
Osxr 038) axa Ox Oax4 O4x4 Osx bey 
Pa + B.(v . v) b? (ptoé, + pt,V,) ~ np, (v : v)v, 
= b! det 
(tp + B.)v" ptob*& + BY — nB,v"v, 
= b[ptgb? — n(v- v)P [Abt + Bb*(v - v) + C(v- v)*] 
= p*rige,b TY (6? —ca(v- 0) (B2) 
a=1,+ 
where, as we have obtained in Eqs. (A4) and (A5), and in 5 4n 
the text below it, B=-t, (veto ar Gar 3 ar ox, — PtptoQ, (B3b) 
C = tp(pc?tg + ox,) — B pet (B3c) 
A= pt,t9, (B3a) = TPIPCSTO 3 3)" 
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n= a = 1, j= [n/(pt;)]. and CL = [(-B ae 
VB? — 4AC)/2A]. It is worth mentioning that the assump- 
tions of Proposition I guarantee that 0 < c,;,c4 < 1. Under 
assumptions (Al), 7 > 0, and conditions (21) in the strict 





Nyt 














form, then one obtains that det(&,A%,) =0 only if 0<c, <1 
bli —Volio 
—vlig blo 
det(E,Aj) = det] —v'lio  Orox10 
—v'Ti9 — O10x:10 
—v1i9 — Ojox10 
(b? — vv, )lio 
—v°I\o 
a aodet -v'lio 
—v"l19 
—v7 119 
= (u%ga)(Sa8")"°. 


Again, note that if & is timelike, then det(é,A%) # 0. This 
completes the proof of (1). 

As for (ii), let us define 6, = €, + Aé, and make the 
changes € > @ in the determinant calculations above. 
Then, the eigenvalues A are obtained from the roots 
of det(p, 2%) = det(p,A%,) det(p,A%)(u%p,)'® = 0. Note 
that the general form of the equations implies that the roots 


ba = —ugu py + Mids obey 


(u“ba)” — BA” babs = 0 (BS) 





Z = P{APC Cp (urg,)? 
> PAP C.Cp(uKE,)° 
= B{(-EEq) (CPEp) ale 


zs APE Es(uNl,) = 
+ APE Es(uhl,) ~ 
[(u%Ea) (uC) 


In the operations above we used the fact that 0 <P < 1, 
(AME nls)” S (AMEE, )(AMC,C,) from the Cauchy- 
Schwarz inequality and that € is timelike and ¢ spacelike. 
Thus, causality guarantees reality of the eigenvalues. 
Now we turn to the problem of completeness of the set of 
eigenvectors. We Deen by counting the apenly ee 
dent eigenvectors of #\”")A%, where pi”) =¢, + AWE, 


and A\”) are the eigenvalues of the matter sector and are 
obtained by means of Eq. (B6) in the cases # = cy = 0 
when a = 0 and f = c, when a = 1, +. Let us define an 
arbitrary vector, 





2(u7Eq) (uP Cg) AMEC, — PU(A% C65) (AXE,E,) > 
2(uE,) (uP) AME, C, ~ 
SAME Ca yO. 
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(with the equality holding only in the case a = 0); ie., the 
equation b2 — c, (vg + vq) = 0 gives €, 4 such that &, ,é% = 
—b? + 04+ Vg = (1—Cq)¥q* Vg > 0. Thus, if ¢ is timelike, 
then (i) is guaranteed for the matter sector as well. As for 
the gravity sector, one obtains that 


—V1\9 —Vlio 3110 
O1ox10 ~~ O10x10 ~— O10x10 
bly ~~ Ot0x10 ~~ O10x:10 
Oiox10 = Lio Ot x10 
Oiox10 10x10 = Lio 
O1ox10 10x10 O10x10 10x10 
blig — Ot0x10  O10x10 ~— O10x10 
Oiox10 = L190 =~ Ot0x10 — O10x10 
Oiox10 O10x10 = PLi0~— O10x:10 
Oiox10 O10x10 O10x10 = BLio 


(B4) 





where, from causality, in any of the above cases we have 
that 0 < # < 1. Then, for each f, the eigenvalues A are 





BAM ECs) — (ua) (u%Sa) & VZ 


A= 
(u%E,)° — BAPEEp 





where, since €&,é% < 0, 
because 0 < # < | and 


then (u%E,)? — BAVE,E, > 0 


(AEC) 1h 


(AMC ap )(AME, 6) + (AME als)" 





ri) _— 


(B8) 
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Then, for each of the eigenvalues A”) a =0,1,+, we must where only 19 out of the 24 components H", J", J 
verify how many of the 29 variables in the vector (B8) are free are free variables because of the 1+ 1+ 3 con- 
parameters under the equation #1") A¢r"") = 0. In fact, this straints Bob H ° : BrPyol *=0, J{=0, and 
is the dimension of the null space of the matrix pa’ f) A® and At agen Jy + AY ger B ne '=0 (mote that the last 
corresponds to the number of linearly independent eigen- four equations are not all independent since the 
vectors of A”). The eigenvectors are the following. contraction with uw, is identically zero, resulting in 

(i) Ay”: This si has muleplicty 19. The eigenvector three independent constraints). Thus, the multiplic- 


that obeys d. m) Agr (m) _ 9 ig ity of Ag equals the number of LI eigenvectors, 
Oa 





ie., 19. 
(ii) A or In this case each of the two eigenvalues have 
0 multiplicity 3 since n; = 3 in Eq. (B2) (note that 
Osx since we assumed here that 7 > 0, then c, # 0 and, 
HE thus, c; # Co and the eigenvalues are different from 
hi the case cy = 0). We may perform some rie 
rr) — oN (B9) row operations over the linear system (”) A¢r\”) = 
0 0 to obtain, by imposing b? — c,(v- v) = 0 @emem- 
Jt ber that b = u*#, and v* = A” pz after the change 
o> p), 
J 





tb? +f,(v-v) bptod, + born, — ev, Og Ora O1xa O1x4 Oia Orn 



































































































































0451 Kv O44 O44 O4x4 04.4 O44 O4x4 
—vh O44 b& Oana Odxa Oana Odxa Oana 
0451 B O44 b O4x4 04.4 O44 O44 A ) = (B10) 
0431 —v9h O44 O4n4 BO) Oayg Oana Oey 
0451 —0 6, Osa O4xa Oana DE Oda Ox 
0451 — v6, Oia Oriya Ota Ona 5bOD Ona 
0451 —030) O44 O4xa Odxa Oaxa Oax4 bey 
| 
where where, from the 29 + 29 = 58 i oea ue = 
above eigenvectors (29 for AM" * and 29 Aj (m) 
K,= ns, + (ver Cra nba) «| cases), they are subjected to the following 26 i 26 
constraints: | + 1 = 2 constraints, 
x [z-b +B,(v , v)| = (tp + B.)[b*pt9é, 
2 _ . m 
+b*pt,v, —np,(v-v)v,]. (B11) [c.b2 + B(vs + v4)\Fs + baptg’ + im \gv 
This enables us to find the eigenvectors, np, (vs - v4) 
+ bypt,veG’ — a a v7 G4, = 0, 
Fy =a 
GY. : + 
. 1+ 1 =2 constraints KG = 0, 4+4=8 con- 
A, straints b.H!. =v'LFs, 4+4=8 constraints 
i rr nv. v+G" + B24 = 0, and the 16 + 16 = is con- 
= = 7 (B12) eae (m) +, (m) 
1 yp |? straints b*J'_, = vf G'., where “gy, = Ay" €4 + 
a ¢, and b* and v% are defined in terms of “go. 
+ Hence, there is a total of 3+ 3 = 6 free parameters. 
i‘. : Once again, the degeneracy equals the number of LI 
“I eigenvectors. 
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(iii) (Ax)*: Since there is no degeneracy in these four 
last eigenvalues and they are distinct from the others 
because cx # 0 in the strict form of the inequalities 
in Eq. (21) and different among them, then one has 
four LI eigenvectors. 

Thus, the system has 19+6+4= 29 LI eigenvectors. 
Therefore, there is a complete set in R**, namely, 


{ri”) }72., such that be Aare” =). Hence, we can use 
the 29 linearly independent set S‘”) = {Ro }29 , to verify 
that 








(B13) 


obeys (€, + Aw”) E,) MeR” =0. 

Now, before we discuss the gravity sector {F'4, F%}, let 
us look at the sector containing the original fields ¢, n, u’, 
and g,,,. In this case, let us define 


(B14) 


4) igs a 161 column vector. Then, (f+ 


where r‘ 
AM E,)UeR = 0 reduces to the eigenvalue problem 
utp ior? = 0 whose eigenvalues are uzg = 0, ie., 
AM = €,u%/E,u*. Thus, the eigenvectors may be any basis 
of R!°. Let {rO}16 | be a basis of R!°. Then, the set S@ = 
{RO} 16 | is a linearly independent set of 16 eigenvectors 
of oa, 

To finalize the eigenvector counting we have to analyze 
the sector containing F, and F%. In this case, let us define 
ee 














O1ox10 O10x10 

WEG Tig (ue G)To 
*p\9 Aa ~ =A'#*g) I O10x10 
A776) T19 —— Otoxt0 
—A3* Pp) T19 —— O1ox10 


which has 40 pivots and 10 independent variables (corre- 
sponding to the variables associated to the first 10 col- 
umns). Thus, there are 10 linearly independent vectors for 





each eigenvalue ~A\: i.e., there is a set Cr, eerie of 





20 linearly independent vectors with corresponding wy, = 
[2 Az }-!L¢*r\ coming from Eq. (Bl6a) such that 


ml 


Si) = RO ~RO}10 | where 
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; (B15) 


O16x1 


where w is some 29 x 1 columns vector while r‘) is a 
50 x 1 columns vector. The eigenvalues of this sector are in 
Eq. (B4) and are given by AD) = u%¢,/wE, coming from 
up = 0 (here #\%), = 6, + AWE,) with multiplicity 30 
and corresponding to # = 0, and the two roots =a) with 


multiplicity 10 coming from pega — 


[upp + Avg “oY, = 0, which corresponds to 


lla 














each 














f= 1, 1e., gravitational waves moving at the speed of 


light. Then, the eigenvalue problem p UeR =0 
reduces to the two equations: 


bDALw, = L&r, (B16a) 


bd 2AGrY = 0. (B16b) 








For the eigenvalues *A\”, one obtains that det| “pag # 
0 because the root / = 1 has been eliminated from the 
matter sector (remember that c, < 1). Thus, there exists a 


solution of Eq. (B16a) for each 79) in Eq. (B16b). One 
needs to count the number of linearly independent rd) for 
AY, i.e., the number of vectors in the basis of the kernel of 
gi. As . In this case, after some elementary row operations 
[look at the second equality in Eq. (B4) after setting 
b? = v- v] one obtains that 























O10x10 10x10 O10x10 
O10x10 10x10 O10x10 
(ub) Tio O10x10 O10x10 (B17) 
O10x10 (ut 619 O10x10 
O10x10 O10x10 (uo I 
Wib 
+ 7 (9) + 
Rhy = a 
O16x1 
is a linearly independent set of 20 eigenvectors 


of nt". 
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As for the eigenvalue A®, note that in this case 


det (9) aa] = 0 because # = cy = 0 is also a root of this 
equation. Thus, for every solution 7) in Eq. (B16b), 
Eq. (B16a) can be either undetermined or have infinite 
solutions. However, for any two different solutions, say, w 

and w2 for one 7), the difference between R‘!)! — R‘9” 
corresponds to a vector in the space spanned by S(”), that 
lies in the kernel of o 9) aA. Therefore, since we are 
counting the number of linearly independent eigenvectors, 
we must choose one particular solution w,, if it exists, for 
each 72), We begin by solving Eq. (B16b). Let {/f = 
u“, 15,15} be a set of linearly independent vectors that are 
Cat AWE, to wit, 76? = 0 and 
{e,}1, be any basis of R'°. Then, one may verify that the 
30 linearly independent vectors 


orthogonal to ge? = 


(B18) 


‘O.ac = c\a 


satisfy dA ay{s) ac — 9. Now we must solve Eq. (B16a), 


where 
013.1 013.1 
oe 9) naa Be), porte 
ole. = vi ore acan = Ky ce Ie ’ (B19) 
Pa os"1(Ca)a ra 
HOB (ex)a Ole 
where we defined 
1 
Ky = 5 ne = Sp) ul uw) (€a)ep 
SS 
Let us look for the particular solution 
0 
Wee = | Om (B20) 
PTOVac 
02051 


Note that 
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O13x1 
(9) bu 
ee 
Bee ye ; 
Bb? he 
BPS A¥ec 


oi! ) A%w ac 


(B21) 


and then, by inserting Eqs. (B19) and (B21) into 


Eq. (B16a), one finds that 


poo ee Kote. (B22) 
This leads to the solution ya. = K,/¢/f, and, thus, 
0 
-K,| 
Wac —= PTO K pL (B23) 
Be ae 
O20x1 
As a consequence, the set So) = = {Rr pr RY : ) RY), — 


Bey Re, Re} with 


Wac 
(9) 
n= | 2, 


O16x1 


is a linearly eee set of 30 ai of G1 


Thus, 6 = S™) US®US (9) U so) contains a Sapte 
set of eigenvectors R of 6,2°R = 0 in R®. This completes 
the proof. : 
We remark that the assumption that the inequalities hold 
in strict form is technical. If equality is allowed, then the 
multiplicity of the eigenvalues might change. This is 
because with equality one can have c, = 0 for a= 1 or 
+ and thus the characteristics defined by b” — c,(v- v) =0 
can degenerate into the characteristics b = 0. Since the 
latter is already present in the system, the multiplicity of the 
characteristics would change. This does not mean that the 
system would not be diagonalizable. Nor does it imply that 
local well posedness, established in the next section, would 
fail [222]. However, a different proof would be needed to 
show diagonalization in the case c, = 0 in the cases a = | 
or +. We believe that treating this very special case here 
would be a distraction from the main points of the paper. 
We also recall that already in the case of an ideal fluid, a 
different approach to local well posedness has to be 
employed when the characteristics degenerate [223]. 
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APPENDIX C: PROOF OF THEOREM II 


As usual in studies of the initial-value problem for 
Einstein’s equations [22], we embed & into R xX and 
work in harmonic coordinates in the neighborhood of a 
point. Observe that we already know the system to be 
causal under our assumptions, thus localization arguments 
are allowed. 

The equations to be studied read 


uXwO gn + nud, 02 gu” +B; (n, u,g)07g = B, (On, Ou, 09), 
(Cla) 


uur wO,OguY + B,(n, e, u, g)0?g = Bo(An, de, Ou, Og), 
(C1b) 


By (uth + AM(@yP))A,Agn + CHPOOge + CPO, gu? 


+ BY(n, e, u, 9)0°g = Bk (On, de, Ou, Ag), (Clc) 
GOO 6 Iu = By,,(On, Oe, Ou, Og), (Cld) 
where 
CY = (00 == ") AM ay) — not 
+ p(t.+ To) ul" A\uP) + toputu’&h, — (C2a) 
GHP — yt (p,A% + 7,u%u) + (B, + tp) AMuP), —(C2b) 


and the notation for the B’s and B’s follow the same 
construction as in Sec. IV. 
We can write Eq. (C1) in matrix form as 


M(I)Y = NIV), (C3) 


where B = (e,n, u’,g,,)' is a 16 x 1 column vector (we 
count only the 10 independent g,,), 8(O'V) is also a 16 x 1 
column vector containing the Jt’s, i.e., the lower-order 
terms in derivatives of each equation, and 


M(@)—-B(@) 


M(O) = ; 
(9) Oi0x6 GF OgOgli0 


(C4) 


The 6 x 10 matrix 6(Q) contains the terms B0?g, while 


0 uruP né(7yP) 
M(0)=]| 0 0 u,utub | O74. (C5) 
Hop Br (ute ote AM(ayP)) cer 
Let us compute the characteristic determinant of the system 


and its roots, ie., det[9t(é)] = det[M(é)|(E%E,)!° = 0, 
where the substitution 0 > & takes place. The pure gravity 
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sector has the roots €*€, = 0. As for the matter sector, by 
again defining b = u%G,, v' = A”’E,, v-v = v'v,, and 


Cf = CO E.es 
= [ropb? — ou: w))8t + (top 6-2) a 


+ p(t. + Tg)bu"v,, (C6a) 


Di = («wp —o- 3 ~ nba) vt, + [topb* — n(v- v)]8, 
(C6b) 


Gt = GEE 


= |B,(v-v) +7,b7|u" + (8. + tp) bv", (C6c) 


where D/ is the same as the one defined in Eq. (A3), we 
obtain that (by carrying out some elementary row operations) 





0 i nbé, 
det[M(é)] =det | 0 0 bu, 
G" B,,[u“(v-v)+bv"| Ct 
b3 
~ Topb?—n(v-v) 
0 b né, 
xdet| 7-b7+£,(v-v) B,(v-v) p(te+t9)bv, |. 
(B.+t+P)bv" B, bv" dD) 


(C7) 


The last determinant is the same as the one obtained in 
Eq. (A4) and the result turns out to be 
det[M(2)] = —b'[prgb? — n(v- 0)? 

x [Ab* + Bb*(v- v) + C(v- v)"] 


= =p tyne.) II (Cates ~ eA 


a=14 





(C8) 
where, as in Eqs. (A5), 


A= pito, (C9a) 


4 
=-t1, (veto +04 5 + ox, —pTptg, (Cb) 





C= Tp(pcztg + ox;) —B, (< + 2) : (C9c) 


while c, = [n/(pt,)] and cz = [(—B + VB? — 4AC)/2A], 


while 7, = 2 and # = 1. Note that the characteristics are 
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still the same as in Sec. IV, as expected, although the 
multiplicity of the roots changed (and there was no reason 
for the multiplicities to be the same). We conclude that the 
characteristic determinant of the system is a product of strictly 
hyperbolic polynomials. We verify at once that the system is a 
Leray-Ohya system [21,224] for which the results of 
Ref. [225] (see also Ref. [105]) apply. Thus, if the initial 
data are quasianalytic (see Ref. [75]), we obtain quasianalytic 
solutions. 

Denote the initial dataset in the theorem by D and let D, 
be a sequence of quasianalytic initial data converging to D 
in H™ (see Ref. [27] for the definition of H). Let W, be 
solutions corresponding to Dy (which exist by the fore- 
going). In order to finish the proof of the theorem, it 
suffices to show that W, has a limit in H’. The limit will 
then be a solution with the desired properties because we 
can pass to the limit in the equations since N > 5. 

According to the arguments given in Sec. 16.2 of 
Ref. [146] or in Refs. [107,108], the diagonalization 
obtained in Sec. VB implies that Y defined in Eq. (24) 
admits a uniform bound in H‘~!, and uniform difference 
bounds in H’~? also holds. We apply these bounds to the 
vector Wy corresponding to ‘Wy. We see at once that the 
uniform H’~! bounds for ¥, imply uniform H” bounds for 
W,, and the difference bounds imply that ‘Wy is a Cauchy 
sequence in H~?, thus converging in this space. But low- 
norm convergence combined with high-norm boundedness 
implies that the limit is in fact in H™ [226]. ] 

We observe that a similar local well-posedness result 
holds for the fluid equations in a fixed background. 

We recall that a standard tensorial argument [22] 
guarantees that the solution established in Theorem II is 
intrinsically defined; i.e., given the data, which are defined 
independently of coordinates or gauge choices, there exists 
a spacetime where Einstein’s equations are satisfied, and 
this spacetime is defined without any reference to coor- 
dinates or gauge choices—even if in the process of proving 
that this spacetime exists one has to work in a specific 
gauge and coordinate system. Therefore, even though we 
used the harmonic gauge in the proof, the existence of the 
solution is guaranteed for other choices as well. This logic 
is similar to showing that a map from a finite-dimensional 
vector space into itself is invertible: one can choose a basis, 
write the matrix of the linear transformation with respect to 
that basis, and compute its determinant. The map is 
invertible if and only if the determinant is nonzero, and 
this conclusion (the invertibility or not of the linear map) is 
independent of any basis choice—even if to show that the 
map is invertible we picked a basis and computed the 
determinant with respect to that basis. 

We note, however, the following subtlety which is very 
relevant for numerical simulations. The fact that a unique 
solution is guaranteed to exist for given initial data, and that 
this solution is well defined regardless of gauge choices, 
does not imply that such a solution can always be 
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reconstructed from an arbitrary gauge. In other words, 
suppose we write the equations in a different gauge. Jf we can 
numerically integrate them, we will obtain the solution found 
in Theorem II written on that gauge (modulo numerical 
accuracy). However, it is possible that the gauge we chose is 
not adequate to solve the equations numerically, so that our 
numerical simulation will not produce a solution. This does 
not mean, of course, that solutions do not exist; it simply 
means that the guaranteed-to-exist solution given by 
Theorem II cannot be accessed from that specific gauge. 
To use again our analogy with determinants, suppose we 
computed the determinant on a basis b, and found it to be 
nonzero, but now we are interested in computing the 
determinant numerically using another basis bp. 
Depending on the basis b, and the numerical algorithm 
we use, this might not be possible, which, of course, does not 
mean that the determinant is zero or ill defined. 

Thus, the practical matter of solving the equations 
numerically is not settled by an abstract existence and 
uniqueness result as Theorem II. Such theorems are 
naturally important as they provide the foundations on 
which numerical investigations can be built; i.e., it makes 
sense to look for solutions numerically because solutions 
do exist. But these theorems do not, in general, point to how 
to recover solutions numerically. That is why there is a 
great deal of work dedicated to writing Einstein’s equations 
in different forms and special gauges, even if basic 
existence results for Einstein’s equations coupled to most 
matter models are known, as reviewed in Refs. [2,29]. 


APPENDIX D: PROOF OF THEOREM III 


From causality one obtains that det(n,A%) # 0 as far as n 
is timelike. Thus, we can rewrite Eq. (35) as 


iQ6Y(K) = —ix(—n,A®)7!C,A%5¥(K) 
— (=n, A*)-'|BSY(K). 














(D1) 


Since the eigenvalue problem (37) contains N linearly 
independent vectors r,, one may write Eq. (37) as 
(—n,A*)"CpA*r, = Agr, (D2) 


and define the N x N invertible matrix R = [r,--- ry] 
whose columns are the eigenvectors 1,...,1¥, and the 
N x N matrix 


t 
L=R!= 
Ty 


where the rows [,, are the left eigenvectors of (—ngA“)E,A’ 
which, consequently, obey [,r, = 6,, (because RL = Jy). 
Then, we can write 
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5¥(K) = RLOY(K) = S°c,(K)tq = Re, (D3) 


where c,(K) = 1,6(K) is a c-number and c is the N x 1 
matrix 


Therefore, Eq. (D1) becomes 





iQRe = —ikRDe — (—n,A*)'BRe, 








(D4) 





where D is the NxWN real diagonal matrix D = 
diag(A,,..., Ay) and, thus, (—n,A%)~'¢,A’R = RD. By 
multiplying Eq. (D4) by c*R™! from the left one obtains that 

iQ\e|? = —ixe’De — e7R7!(—n,A*)"'BRe. (D5) 


Since D is real and diagonal (which gives c’Dc € R), 
Q=y,(-iT + c’k)), and x = y¢(-id/cT + d/k;), then 














Dre! (y,ly + yed'c;D)e = —Re{ctR-!(n,A%)~'BRe]. (D6) 


On the other hand, note that y,J,y + yed! c;D is diagonal 
with elements 


(Yn + yed!cjD) aq =Yn cs yedicjNg >0 (D7) 


because |d/c;| <|c'| <1, A€[-1, 1], and y, > y¢ from 
Eq. (36). Hence, y,Jy + red! c;D is a positive Hermitian 
matrix and ¢'(y,Jy + red! c;D)c > 0. The consequence is 
that Cp < 0 if and only if 














Re[c'R7!(n,A*)'!BRe] > 0. (D8) 

Now, let O be the LRF and ©’ some other boosted frame. 
The connection between the two frames is given by the 
Lorentz transform ¢/ = y(t — v'x;), x = v(x —v't), and 
x'!, = x!, where || and L stand for the components parallel 
and perpendicular to v', respectively. This can be com- 
pactly written as X’“ = AY X”. Thus, one obtains that K’" = 
ALK’ and &¥’(K') = M6¥(K) from the structure of 
Eq. (33) (where M is an N x N invertible matrix), leading 
to A’ = AYMA“M~! and B’ = MBM~!. In particular, 
CgA*% = M'CLA’*M and nzA® = M7!(n,A’*%)M. From 
Eq. (37), these relations give R’ = MR, with the same 
eigenvalue A in both frames. Then, since 6¥(K) = Re and 
OY'(K’) = R'c’ = MR, one concludes that c=’, ice., 
cl,(K") = cq(K). Therefore, one arrives at the following 
identity: 
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CTR’ (=n, A'2)“! B’R’¢! _ ciR7(—n, A%)"! BRe. 














(D9) 








However, if the system is stable in the LRF, then Eq. (D8) 
holds and, from Eq. (D9), one automatically obtains that 
Ip < 0, proving that the system is also stable in any other 
frame O’ obtained via a Lorentz transformation. . 
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